solving a system of 4 second order differential equations

i have the following first order equations

z1=x, z2=dx/dt, z3=y and z4=dy/dt

which relate to the following second order equations

x'=z2,x''=4*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3), y'=z4 and y''=2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)]

i have created the following function file to slove the system using ode45. however i get the error index exceeds number of array elements.

function zprime=secode(t,z)

%secode: Computes the derivatives involved in solving the first order

%system

zprime=[z(2);2*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3);z(4);-2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)];

this is my solution file code

xspan=input('please enter your time range:');

initp1=input('please enter initial condition for x(0)=?:');

initp2=input('please enter initial condition for y(0)=?:');

initv1=input('please enter the initial condition for dx/dt(0)=?:');

initv2=input('please enter the initial condition for dy/dt(0)=?:');

cond1=[initp1 initp2];

cond2=[initv1 initv2];

[x y]=ode45(@secode,xspan,cond1,cond2);

re=sqrt((z(1)+E)^2+z(3)^2);

rm=sqrt((z(1)-M)^2+z(3)^2);

M=1-E;

E=0.0123;

any help would be greatly appriciated.

[darova]

Can you atach original equations? Paper or something

[darova]

Can you show also this?

[darova]

Are you running out of space?

Better be to write as rows

[Stephen23]

Original question from Google Cache (and formatted properly):

"solving a system of 4 second order differential equations"

i have the following first order equations

z1=x, z2=dx/dt, z3=y and z4=dy/dt

which relate to the following second order equations

x'=z2,x''=4*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3), y'=z4 and y''=2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)]

i have created the following function file to slove the system using ode45. however i get the error index exceeds number of array elements.

function zprime=secode(t,z)
%secode: Computes the derivatives involved in solving the first order
%system
zprime=[z(2);2*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3);z(4);-2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)];
this is my solution file code
xspan=input('please enter your time range:');
initp1=input('please enter initial condition for x(0)=?:');
initp2=input('please enter initial condition for y(0)=?:');
initv1=input('please enter the initial condition for dx/dt(0)=?:');
initv2=input('please enter the initial condition for dy/dt(0)=?:');
cond1=[initp1 initp2];
cond2=[initv1 initv2];
[x y]=ode45(@secode,xspan,cond1,cond2);
re=sqrt((z(1)+E)^2+z(3)^2);
rm=sqrt((z(1)-M)^2+z(3)^2);
M=1-E;
E=0.0123;

any help would be greatly appriciated.

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[darova]

I believe that re and rm should be functions too

re = @(z) sqrt((z(1)+E)^2+z(3)^2);
rm = @(z) sqrt((z(1)-M)^2+z(3)^2);
zprime = @(t,z) [z(2)
    2*z(4)+z(1)-((M*(z(1)+E))/re(z)^3)-((E*(z(1)-M))/rm(z)^3)
    z(4)
    -2*z(2)+z(3)-((M*z(3))/re(z)^3)-((E*z(3))/rm(z)^3)];
tspan = [0 5e-3];
z0 = [0 1 0 1];
[t, z] = ode45(zprime,tspan,z0);

Pay attention

Answer 2

[Steven Lord]

The third input argument to ode45 is the vector of initial conditions and the fourth is the options. You're trying to pass vectors of initial conditions in as the third and fourth inputs, and by looking at the third input ode45 thinks you only have two ODEs not four. That means it's going to call your ODE function with a two-element vector z. Combine the initial condition vectors into one vector and pass the combined vector as the third input.