MATLAB: Doubt in MATLAB coding.

2nd orderdifferential equations

Hello. I'm trying to plot a differntial equation using MATLAB. When I run my program, I'm getting different graph.
I've attached the differential equation and my MATLAB coding below.
Can anyone tell me how to use the boundary conditions correctly?
function abcd
ex4init=bvpinit(linspace(0.00001,100,21),[0,0]);
sol=bvp4c(@ex4ode,@ex4bc,ex4init)
plot(sol.x,sol.y(1,:),'blue');
end
function f=ex4ode(r,y)
a=0.001391304348;
b=19.31459475;
f=[y(2)
    -(2/r)*y(2)-a*(1-exp(b*y(1)))
        ];
end
function res=ex4bc(ya,yb)
res=[
    ya(1)+1
     yb(2)
            ]
end
This is the graph which I got.
But this is the exact graph.

Best Answer

Hi Joy Salomi,
Thanks for the paper. The differential equation should be written in spherical coordinates (edit: I think there is a mistake to assume R to be constant). Since R is assumed to be constant there is no singularity. I adjusted the code and tested it in Matlab 2019b. Depending on chosen xmesh results might differ. You have to try a bit to get a stable solution. I attached an example for the D = 200 nm graph.
abcd.m
function abcd(ah,D,phi_s)
ex4init=bvpinit(linspace(0,D/2,3e3),[phi_s 0]);
sol=bvp4c(@ex4ode,@ex4bc,ex4init);
plot(ah,sol.x*1e9,sol.y(1,:));
ah.XLabel.String = 'r [nm]';
ah.YLabel.String = 'phi [V]';
    
    function dphidr=ex4ode(r,phi)
        a = -1.602176634e-19*1e17*(1e-2)^(-3)/(8.8541878125e-12*11.5);
        b = 1.602176634e-19/(physconst('Boltzmann')*600);
        dphidr = [ phi(2)
                   a * ( 1-exp( b*phi(1) ) ) ];
    end
    function res=ex4bc(ya,yb)
    res=[ ya(1)+ phi_s
          yb(2) ];
    end
end
Example call:
fh = figure;
ah = axes(fh);
hold(ah,'on')
arrayfun(@(dIn) abcd(ah,200e-9,dIn),[0.1,0.3,0.7,1])
Result:
Kind regards,
Robert
Related Question