1. The problem statement, all variables and given/known data There is a grounded sphere of radius R in the origin of the coordinate system. In the distance L (L>R) from the sphere’s center there is a point charge Q. The electric field (both intensity and potential) should be computed in the area of radius rg = 5L (in the plane containing the sphere’s center and the charge). The solution should be delivered in Scilab. 2. Relevant equations Poisson and Laplace equations solved numerically. ∇2 [itex]\varphi[/itex] = 0 Laplace ∇2 [itex]\varphi[/itex] = f Poisson 3. The attempt at a solution Using the Taylor series expansion of derivatives, one can convert the differential equation to the difference equation. In the case of the region outside the sphere, where the Laplace equation is used, basically, the electric potential in a given point is computed as the arithmetic average from its four neighboring potentials, the north, west, east and south potential. However, this is not true for the region belonging to the point charge, where I assume, the Poisson equation should be used. Can you help me finding a way to use the Poisson equation in the code? Can you take a look at the code to check if until know everything is correct? Are my assumptions with Laplace and Poisson equations right? My code so far: clear // clear all the variables present in the memory xdel(winsid()) //close all the graphical windows multiplier = 3; //used to adjust the density of the grid density = 30 * multiplier; //grid's density L = 8; Rg = 2*L; Rc = 0.5; R = 5; q = 10; X0 = - Rg; Xk = Rg; Y0 = - Rg; Yk = Rg; x = linspace(X0, Xk, density); y = linspace(Y0, Yk, density); [X,Y] = ndgrid(x,y); dx = x(2) - x(1); dy = y(2) - y(1); sphere = zeros(density, density); r = sqrt(X.^2 + Y.^2); charge = sqrt ((X-L).^2 + Y.^2); sphere(find(r <= R)) = 0; sphere(find(charge <= Rc)) = q; //xset('colormap', jetcolormap(256)) //plot3d1(x,y,sphere) dok=1e-3; count=0; jeszcze=%t; while jeszcze do count = count + 1; nu=([sphere(:,$), sphere(:,1:$-1)]+[sphere(:,2:$),sphere(:,1)]+[sphere(2:$,:);sphere(1,:)]+[sphere($,:);sphere(1:$-1,:)])/4; nu(find(r <= R)) = 0; nu(find(charge <= Rc)) = q; nu(:,1)= (4*nu(:,2)-nu(:,3))/3; nu(:,$)= (4*nu(:,$-1)-nu(:,$-2))/3; nu(1,:)= (4*nu(2,:)-nu(3,/3; nu($,:)= (4*nu($-1,:)-nu($-2,/3; delta = max(abs(sphere-nu)); sphere = nu; if (delta < dok) then jeszcze = %f; end; end; xset('colormap', jetcolormap(64)) plot3d1(x,y,sphere) xtitle(msprintf('%g loops were performed, max incerement = %g / %g',count,delta,dok),'x','y') grx = -(sphere(3:$,2:$-1) - sphere(1:$-2,2:$-1))/2/dx; gry = -(sphere(2:$-1,3:$) - sphere(2:$-1,1:$-2))/2/dy; dlgr = sqrt (grx.^2+gry.^2); maxdl=10; wspdl=dlgr/maxdl; wspdl(wspdl<1) = 1; ngrx = grx./wspdl; ngry = gry./wspdl; scf(); contour2d(x,y,sphere,10) xtitle("","x","y"); //scf(); co = 2; champ(x(2:co:$-1), y(2:co:$-1), ngrx(1:co:$,1:co:$), ngry(1:co:$,1:co:$)); a =gca(); a.isoview = 'on'; // check the correctness of the solution ip=24*multiplier; ik=26*multiplier; jp=15*multiplier; jk=22*multiplier; plot2d(x([ip,ip,ik,ik, ip]), y([jp,jk,jk, jp, jp]), style = 5) gr1 = -(sphere(ip:ik,jk+1) - sphere(ip:ik,jk-1))/2; gr1(1) = gr1(1)/2; gr1($)=gr1($)/2; gr2 = -(sphere(ik+1, jp:jk) - sphere(ik-1,jp:jk))/2; gr2(1) = gr2(1)/2; gr2($)=gr2($)/2; gr3 = (sphere(ip:ik,jp+1) - sphere(ip:ik,jp-1))/2; gr3(1) = gr3(1)/2; gr3($) = gr3($)/2; gr4 = (sphere(ip+1,jp:jk)-sphere(ip-1,jp:jk))/2; gr4(1)=gr4(1)/2; gr4($)=gr4($)/2; F1 = sum(gr1); F2 = sum(gr2); F3 = sum(gr3); F4 = sum(gr4); FLUX = F1+F2+F3+F4; xtitle(msprintf('Number of divisions is %g x %g Net flux = %g \nF1=%g F2=%g F3=%g F4=%g',density,density,FLUX,F1,F2,F3,F4),'x','y') You can try it out on your own.