if you read the help on ode45, it is pretty clear on how to do this.
The equations you wrote are 4 separate pairs of second order equations that do not affect each other.
I don't understand the y function on the right hand side. You show it as a vector of 8 elements. y(t) should be a continuous function of time. Please clarify.
Since the 4 systems operate independently, you can solve this with one call to ode45, using a solution vector with 8 elements, or you can do it with 4 calls to ode45(), with a solution vector of two elements each time.
Doing it with one call, you will have the solution vector:
x(1..4)=x1,x2,x3,x4. x(5..8)=x1',x2',x3',x4'
Write a function that returns the derivatives of the 8 elements. Pass that function to ode45. You save this function as a separate .m file, or you include it after your main code.
For example
function dxdt=odefun(t,x)
dxdt=zeros(8,1); %force x to be a column vector
dxdt(1)=x(5);
dxdt(2)=x(6);
dxdt(3)=x(7);
dxdt(4)=x(8);
dxdt(5)=-3.797*x(5)-1442*x(1);
dxdt(6)=-10.815*x(6)-11698*x(2);
dxdt(7)=-16.1997*x(7)-26227*x(3);
dxdt(8)=-19.1621*x(8)-36719*x(4);
end
In the code above, I have not include the y(t) term since I do not understand your explanation of it.
You will need to specifiy initial conditions for the 8 elements of x() in your main code..
