Code not giving correct result when constants are given of the order of 10^-10, whilst when i give constant of the order of 10, it runs.

Following is the code for solving coupled diff eqn using rk4, it is for pupulation dynamics of a two level laser system. Here the real values of h=0.02e-9, xn=2e-9, A=5e-15, B=2e-9. With these the code does not furnish proper result, but with the following valued, it runs.

Also i don't know how to use vector function handle hence the lengthy version.

clc
clf
syms h xn A B ;
Inputh=(0.02);
Inputxn=(2);
InputA=(5);
InputB=(2);
h=vpa(Inputh)

h = 

0.02

xi=0;
xn=vpa(Inputxn);
x=vpa(linspace(xi,xn,100));
A=vpa(InputA);
B=vpa(InputB);
f_1=@(y,z)(-A*y+A*z+z/B)

f_1 = function_handle with value:
    @(y,z)(-A*y+A*z+z/B)

f_2=@(y,z)(A*y-A*z-z/B)

f_2 = function_handle with value:
    @(y,z)(A*y-A*z-z/B)

y=zeros(1,length(x));
z=zeros(1,length(x));
y(1)=100;
z(1)=0;

for i=1:(length(x)-1)
    k_1=h*(f_1(y(i),z(i)));
    l_1=h*(f_2(y(i),z(i)));
    
    k_2=h*(f_1(y(i)+k_1/2,z(i)+l_1/2));  
    l_2=h*(f_2(y(i)+k_1/2,z(i)+l_1/2));
    
    k_3=h*(f_1(y(i)+k_2/2,z(i)+l_2/2));
    l_3=h*(f_2(y(i)+k_2/2,z(i)+l_2/2));
    
    k_4=h*(f_1(y(i)+k_3,z(i)+l_3));
    l_4=h*(f_2(y(i)+k_3,z(i)+l_3));
    
    K=(k_1+2.*k_2+2.*k_3+k_4)./6;
    L=(l_1+2.*l_2+2.*l_3+l_4)./6;
    
    y(i+1)=y(i)+K;
    z(i+1)=z(i)+L;
end
hold on
plot(x,y)

plot(x,z)
hold off

[vipul kumar]

Please suggest why is it not giving the same graph when i pass the parateres of the order of 10^(-9).

Answer 1

[Sandeep]

Hi Vipul kumar,

It is my understanding that you could obtain the Graph given by using the values without exponential notation but the same graph is not produced when the values used are in the order of 10^(-9).

This discrepancy is due to the fact that both the graphs use the same X-limit and Y-limit and points plotted in the order of 10^(-9) become insignificant when the Y-axis limit is in the order of 10. You can refer the following documentation page to get more insight on setting custom X and Y limits.

Set or query x-axis limits - MATLAB xlim (mathworks.com)

Set or query y-axis limits - MATLAB ylim (mathworks.com)