i want to solve the four coupled Diff Eqs given in code and evaluate it up to Tmax. In last i want to plot rho against z. I expanded the variables with taylor series and apply Runge kutta second order, but i could not get the required plot.

%%%%calculate the trajectory of electron in the magnetic field

%%My Equations

% drho/dt=P_rho (1st Eq)

% dZ/dt=P_z (2nd Eq)

% dP_rho/dt=-rho/(rho^2+Z^2)^3/2 -2*diama*rho (3rd Eq)

%dP_z/dt=-Z/(rho^2+Z^2)^3/2 (4th Eq)

%%%%%%define constant%%%

B = 5.96; %(T)

B_au = 5.96/(2.35*10^5); % in a.u.

Tc = 2*pi/B_au ; %a period0

Tmax = 10*Tc; %the max of trajectory time

diama = B_au^2 /8 ; %the coefficient of diamagnetic term

n = 10^8; %count number

dt = Tmax / n; %the time step

rb = 50; %in u.a. (a0)

%%%%%%%define variable%%%%%

z = zeros(1, n +1);

rho = zeros(1,n +1); %initial z and rho

p_z = zeros(1, n+1);

p_rho = zeros(1, n+1); %intial momentum

%%%%%%%%algorithm for calculate z and rho%%%%%%

for

j = 41.5: 0.01:41.6

theta = j/90 *pi/2 ; %intial angle should vary between (0 ,pi/2)

p_r = (2/ rb- 1/4 *B_au^2*rb^2*sin(theta)^2)^(1/2); %intial momentum in r

coordinate

z(1) = rb *cos(theta);

rho(1) = rb* sin(theta);

p_z(1) = p_r*cos(theta);

p_rho(1) = p_r*sin(theta); %%initial conditions

i = 1;

 while 
    i <n+1  && rho(i)^2 + z(i)^2 >= rb^2

    rho(i+1) = 
rho(i)+ dt*p_rho(i)+1/2(-1*rho(i)/(rho(i)^2+z(i)^2)^(3/2)-2*diama*rho(i))*dt^2;

    z(i+1) = z(i) +dt*p_z(i)+1/2*(-1*z(i)/(rho(i)^2+z(i)^2)^(3/2))*dt^2;

    %%using 2nd Runge-Kutta methods to calculate trajectories%%%%

    %%%next step to calculate the coefficients of 2nd term.%%%

    p_z_tt = -p_z(i)*(rho(i)^2+z(i)^2)^(-3/2)+3/2*z(i)*(rho(i)^2+z(i)^2)^(-5/2)*
(2*rho(i)*p_rho(i)+2*z(i)*p_z(i));

    p_z_tz = -(rho(i)^2+z(i)^2)^(-3/2)+3/2*2*z(i)^2*(rho(i)^2+z(i)^2)^(-3/2);

    p_rho_tt =  -p_rho(i)*(rho(i)^2+z(i)^2)^(-3/2)-2*diama*p_rho(i)+(3/2)* rho(i)*(rho(i)^2+z(i)^2)^(-5/2)*(2*rho(i)*p_rho(i)+2*z(i)*p_z(i));

    p_rho_trho = -(rho(i)^2+z(i)^2)^-3/2+3/2*2*rho(i)^2*(rho(i)^2+z(i)^2)^(-3/2)-2*diama;

    p_rho_t = -1*rho(i)/(rho(i)^2+z(i)^2)^(3/2)-2*diama*rho(i);

    p_rho(i+1) =p_rho_t *dt+p_rho(i)+1/2*dt^2*(p_rho_tt+p_rho_trho*p_rho_t);

    p_z(i+1) = (-1*z(i)/(rho(i)^2+z(i)^2)^(3/2))*dt+p_z(i)+1/2*dt^2*(p_z_tt+p_z_tz*(-1*z(i)/(rho(i)^2+z(i)^2)^(3/2)));

    i = i +1;

 end

% figure(1);

% plot(z,rho);% hold on ;

if i ~= n+1 && abs(p_rho(i)*z(i)/rho(i)-p_z(i))<=10^(-3)

      figure(1);

      plot(z,rho);% hold on ;

%print (1,'-deps ', num2str (j*10000));

      i,j

end

end

[John D'Errico]

PLEASE learn to format your code to be readable As it is, you can see only a strung out, unreadable mess here. So read this:

https://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup#answer_18099

Next, exactly what is your question? What do you need?

If you want someone to help you, then make it easy for them to be of help.