Runge Kutta 4 order coupled differential equation

we need the values of z, r, m11 and m1 by solving the above equations but not getting the correct values. Please help me find the error in the code. We are getting repeated value of r, z, m1 and m11.

a = eos1 ;

N = a(:,1) ; % w(:,1) number density

ed = a(:,2) ;% w(:,2) energy density

p11 = a(:,3) ;%w(:,3) parallel pressure

p1 =a(:,4) ; % w(:,4) perpendicular pressure

n = 14600;

r(1) = 10^-28 ;

z(1) = 10^-28 ;

m11(1)= 10^-28;

m1(1)= 10^-28;

f1 = @(z,r,ed,p11,m11) -(r.^2+z.^2).^7 *(4*pi*p11*r) ./(2*z*m11.^3 *r) ./(p11*(r+1)-r*(ed+p11)*(r.^2+z.^2).^6); %dz/dp11

f2 = @(z,r,p1,ed,m1) -m1.^3 *r.^2 ./(ed+p1) ./(r.^2+z.^2).^6 ./(4*pi*((ed+p1)-p1*(r-m1))) ; % dr/dp1

f3 = @(r,ed) 4*pi*r.^2 *ed ./3 ; % dm11/dz

f4 = @(r,ed,z) 8*pi*r*ed*z ./3 ; % dm1/dr

for i = 1:n-1

h = p11(i+1,:) - p11(i,:) ; %step size for parallel pressure

j = p1 (i+1,:) - p1(i,:); %step size for perpendicular pressure

if i ==1

k = z(1) - 0 ; % step size for z component

l = r(1) - 0 ; %step size for r component

else

k = z(i) -z(i-1) ;

l = r(i) - r(i-1) ;

end

w1 = f1(p11(i,:),z(i),r(i),m11(i),ed(i)) ; %%%% for first equation

w2 = f1(p11(i,:)+0.5*h,z(i)+0.5*w1,r(i)+0.5*w1,m11(i)+0.5*w1,ed(i)) ;

w3 = f1(p11(i,:)+0.5*h,z(i)+0.5*w2,r(i)+0.5*w2,m11(i)+0.5*w2,ed(i)) ;

w4 = f1(p11(i,:)+h,z(i)+w3,r(i)+w3,m11(i)+w3,ed(i)) ;

x1 = f2(p1(i,:),z(i),r(i),m1(i),ed(i)) ; %%%%%%%% for 2nd equation

x2 = f2(p1(i,:)+0.5*j,z(i)+0.5*x1,r(i)+0.5*x1,m1(i)+0.5*x1,ed(i)) ;

x3 = f2(p1(i,:)+0.5*j,z(i)+0.5*x2,r(i)+0.5*x2,m1(i)+0.5*x2,ed(i)) ;

x4 = f2(p1(i,:)+j,z(i)+x3,r(i)+x3,m1(i)+x3,ed(i)) ;

y1 = f3(ed(i),r(i)) ; %%%%%%%%% for 3rd equation

y2 = f3(ed(i),r(i)+0.5*y1) ;

y3 = f3(ed(i),r(i)+0.5*y2) ;

y4 = f3(ed(i),r(i)+y3) ;

v1 = f4(r(i),ed(i),z(i)) ; %%%%%%%% for 4th equation

v2 = f4(r(i)+0.5*l,ed(i),z(i)+0.5*v1) ;

v3 = f4(r(i)+0.5*l,ed(i),z(i)+0.5*v2) ;

v4 = f4(r(i)+l,ed(i),z(i)+v3) ;

z(i+1) = z(i) + (1/6)*h*(w1+2*w2+2*w3+w4);

r(i+1) =r(i) + (1/6)*j*(x1+2*x2+2*x3+x4) ;

m11(i+1) = m11(i) +(1/6)*k*(y1+2*y2+2*y3+y4);

m1(i+1) = m1(i) + (1/6)*l*(v1+2*v2+2*v3+v4) ;

end

[Torsten]

Didn't you read my response to your question ? Your update Runge-Kutta formulae are wrong.

They would only be correct in the case that the first differential equation only depends on the first unknown, the second differential equation only depends on the second unknown and so forth. This is not the case for your system.

[Kishan Bajpai]

So sir , How can we apply Range kutta 4th order to solve this type of equation? Is there any way ? @torsten sir

[Torsten]

Let us first analyze your system of equations.

It looks as if your equations (1) and (2) can directly be integrated to give

M_|| (r,z) - M_|| (r,z0) = 4/3*pi*r^2*eps_T * (z-z0)
M_|_ (r,z) - M_|_ (r0,z) = 4/3*pi*z*eps_T* (r^2-r0^2)

Is that correct ?

[Kishan Bajpai]

Actually No , because There is problem that there are all unknown like mass (perpendicular and parallel) and r , z .so How can we apply range kutta 4th order . We have only one condition there r(0) z(0) m1(0) m11(0) are equal to zero

[Torsten]

r and z are the independent variables and thus are not unknown.

Or does your system come from characteristic equations for partial differential equations for P_|| and P_|_ ?

Then please include this original PDE system.

[Kishan Bajpai]

Yes we have to calculate r and z from there , in that equation p11 and p1 are known . Then We have to use that data for m11 and m1 . For p11 ,p1 epsilon i have data of 14601×1

[Kishan Bajpai]

Whatever I wrote the code , when I run it , I get this data .but the problem was there were all same value of r and z and same for m1 and m11

[Torsten]

I can't give advice since from what you write, I don't understand what you are trying to do.

The only thing I can tell you is that your Runge-Kutta formulae are wrong.

Look at

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

to see how to update correctly.