How to simulate using ode23 with multidimentional second order differential equations

!!!At the very bottom, there is a much simpler version of this question!!!

I am trying to model a beam using FEM in Matlab. I understand how to use ode23 for single dimentional problems however, my problem is 10 dimentional and second order. It looks like [M]*(xdotdot)+[K]*x=Q where M and K are 10x10 matrices. I think I have my code close to being right but I'm not quite sure what is wrong. When I debug the code it stops when It gets to the ode function at the bottom. Note that it is calling function sys.m

Here is my code:

clc; clear
format shortE
t=1:.01:10; % time
% -------------------------------------------------------
%{
% User input
l = 3; % length in inches
EI1 = 5*10^6; % EI
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = 200; % spring constant 1
k2 = k1; % spring constant 2
m = 268.2; % mass per unit length
po = 100; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = 75; % point load
%}
% -------------------------------------------------------

% -------------------------------------------------------
% User input
% Validation code
l = 3; % length in inches
EI1 = 5*10^6; % EI
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = 0; % spring constant 1
k2 = k1; % spring constant 2
m = 268.2; % mass per unit length
po = 100; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = [75]'; % point load
zeta = .5;
% -------------------------------------------------------

% -------------------------------------------------------
% Phi function handles
phi1 = @(zeta) 1-3*zeta^2+2*zeta^3;
phi2 = @(zeta) l*zeta-2*l*zeta^2+l*zeta^3;
phi3 = @(zeta) 3*zeta^2-2*zeta^3;
phi4 = @(zeta) -l*zeta^2+l*zeta^3;
Phi = @(zeta)[...
    phi1(zeta)*phi1(zeta) phi1(zeta)*phi2(zeta) phi1(zeta)*phi3(zeta) phi1(zeta)*phi4(zeta);...
    phi2(zeta)*phi1(zeta) phi2(zeta)*phi2(zeta) phi2(zeta)*phi3(zeta) phi2(zeta)*phi4(zeta);...
    phi3(zeta)*phi1(zeta) phi3(zeta)*phi2(zeta) phi3(zeta)*phi3(zeta) phi3(zeta)*phi4(zeta);...
    phi4(zeta)*phi1(zeta) phi4(zeta)*phi2(zeta) phi4(zeta)*phi3(zeta) phi4(zeta)*phi4(zeta) ];

%PPhi = Phi(zeta)
%}
% This function needs to be moved to where it is needed
% -------------------------------------------------------

% -------------------------------------------------------
% Displacement function handles
Qpo = @(po,l) [po*l/2 po*l^2/12 po*l/2 -po*l^2/12]';
QP = @(zeta,P) [-P*phi1(zeta) -P*phi2(zeta) -P*phi3(zeta) -P*phi4(zeta)]';
Qspring3 = @(k1,zeta) k1*Phi(zeta);
Qspring4 = @(k1,zeta) k1*Phi(zeta);
% This function needs to be moved to where it is needed
% -------------------------------------------------------

% -------------------------------------------------------
% Calculation simplification
a=zeros(4,6);
aa=zeros(4);
aaa=zeros(4,2);
I4=eye(4);
A1=[I4 a];
A2=[aaa I4 aa];
A3=[aa I4 aaa];
A4=[a I4];
% -------------------------------------------------------

% -------------------------------------------------------
% load transformations
% p



% P
% -------------------------------------------------------


% -------------------------------------------------------
% Mass matrix, Stiffness matrix
M1 = (m*l/420)*[...
    156 22*l 54 -13*l;
    22*l 4*l^2 13*l -3*l^2
    54 13*l 156 -22*l
    -13*l -3*l^2 -22*l 4*l^2];
M2 = M1; % note, this can be different
M3 = M1; % note, this can be different
M4 = M1; % note, this can be different

K1 = (EI1/l^3)*[...
    12 6*l -12 6*l
    6*l 4*l^2 -6*l 2*l^2
    -12 -6*l 12 -6*l
    6*l 2*l^2 -6*l 4*l^2];
K2 = (EI2/l^3)*[...
    12 6*l -12 6*l
    6*l 4*l^2 -6*l 2*l^2
    -12 -6*l 12 -6*l
    6*l 2*l^2 -6*l 4*l^2];
K3 = (EI3/l^3)*[...
    12 6*l -12 6*l
    6*l 4*l^2 -6*l 2*l^2
    -12 -6*l 12 -6*l
    6*l 2*l^2 -6*l 4*l^2];
K4 = (EI4/l^3)*[...
    12 6*l -12 6*l
    6*l 4*l^2 -6*l 2*l^2
    -12 -6*l 12 -6*l
    6*l 2*l^2 -6*l 4*l^2];
% -------------------------------------------------------
% calculations
K3 = K3 + Qspring3(k1,l/d2); % Khat
K4 = K4 + Qspring4(k1,l/d3);  % Khat
Qpo4 = Qpo(po,l);
QP4 = QP(l/d1,P);
syms Q11 Q21 Q31 Q41
syms Q12 Q22 Q32 Q42
syms Q13 Q23 Q33 Q43
syms Q14 Q24 Q34 Q44
Q1 = [Q11 Q21 Q31 Q41]';
Q2 = [-Q31 -Q41 Q32 Q42]';
Q3 = [-Q32 -Q42 Q33 Q43]';
Q4 = [-Q33 -Q43 0 0]'+QP4+Qpo4;
Q11 = A1'*Q1;
Q22 = A2'*Q2;
Q33 = A3'*Q3;
Q44 = A4'*Q4;
Q = Q11+Q22+Q33+Q44;
Q(1) = 0; Q(2) = 0;
% -------------------------------------------------------

% -------------------------------------------------------
% global calculations
M = A1'*M1*A1+A2'*M2*A2+A3'*M3*A3+A4'*M4*A4;
K = A1'*K1*A1+A2'*K2*A2+A3'*K3*A3+A4'*K4*A4;
Minv = inv(M);
% -------------------------------------------------------

%qstatic = inv(K)*Q

time = (0:.001:22.5)';
x0 = [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]; %[position,velocity]
[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt]=ode23(@(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt) sys(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt,Minv,K), time, x0);

Then my sys.m looks like:

function [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt] = sys(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt,Minv,K)
f = zeros(20,1);
f(1) = x(2);
f(2) = x(3);
f(3) = x(4);
f(4) = x(5);
f(5) = x(6);
f(6) = x(7);
f(7) = x(8);
f(8) = x(9);
f(9) = x(10);
f(10) = x(11);


 delta = -Minv*K*[x(12) x(13) x(14) x(15) x(16) x(17) x(18) x(19) x(20) x(21)]'%+Minv*Q;
 f(11) = delta(1)
f(12) = delta(2)
f(13) = delta(3)
f(14) = delta(4)
f(15) = delta(5)
f(16) = delta(6)
f(17) = delta(7)
f(18) = delta(8)
f(19) = delta(9)
f(20) = delta(10)
end

Here is the easier version of the code:

clc; clear
figure
time = (0:.001:22.5)';
M = [1 2; 3 4];
K = [5 6; 7 8];
x0 = [0;0;0;0]; %[position,velocity]
[t1,x1,x2,x3,x4]=ode23(@(t1,x) trick(t1,x,M,K), time, x0);

plot(t1,x(:,1));
xlabel('Time(sec)');
ylabel('Displacement');
title('Stepped Response(Underdamped)');

Here is the code that it calls:

%2x2 sys
function f = trick(t,x,M,K)
f = zeros(4,1);
f(1) = x(2);
f(2) = x(3);
MinvnegK = -inv(M)*K;
delta = MinvnegK*[x(4) x(5)]';
f(3) = delta(1);
f(4) = delta(2)

[madhan ravi]

Sorry had to close the other question because this question is the same as the other one .

Answer 1

[madhan ravi]

replying to your comment:

clc; clear
figure
time = (0:.001:22.5)';
M = [1 2; 3 4];
K = [5 6; 7 8];
x0 = [1;0;0;0]; %[position,velocity]
[t1,x]=ode23(@(t,x) trick(t,x,M,K), time, x0);
plot(t1,x(:,1));
xlabel('Time(sec)');
ylabel('Displacement');
title('Stepped Response(Underdamped)');

%2x2 sys
function f = trick(t,x,M,K)
f = zeros(4,1);
f(1) = x(1);
f(2) = x(2);
MinvnegK = -inv(M)*K;
delta = MinvnegK*[x(3) x(4)]';
f(3) = delta(1);
f(4) = delta(2);
end

[Jesse Crotts]

Thanks for the answer. Do you know of a faster way to write the values at the end of the function. I know its easy for a 2x2 matrix but if you had a 100x100 matrix what would be the best way to do it:

delta = MinvnegK*[x(3) x(4)]';
f(3) = delta(1);
f(4) = delta(2);

[madhan ravi]

delta(:,100) just use indexing to store the values

Answer 2

[Bob]

It is late :( but might be useful :)

clear; clc;         % clear the desk
t = [ 0 , 10 ] ;    % time interval, use the default interior points
x = [1 0, 0 0]';    % starting point [position, velocity]
B = inv([1 2; 3 4]) * [5 6; 7 8] ;	% define matrix B = -inv(M) * K
s = ode23 (@(tau,z) [z(1:2); B*z(3:4)], t, x);	% solve with ode23()
plot(s.x, s.y(1,:));% check the shape of the solution and its boundaries

disp(sprintf("\tt in [ %d, %d ] \ts in [ %g, %g ]",t,min(s.y(1,:)), max(s.y(1,:))));

	t in [ 0, 10 ] 	s in [ 1, 21846.1 ]