MATLAB: Solving second order ode problem with ode 45

ode45second order ode

I am analysing a mass spring damper system too, but mine has multiple degrees of freedom. I have a function that creates a column vector of dydt once time and a state vector are provided and it works fine. if I use it as funcion in ode45 it has problems. Which is the correct form of input vectorial function for ode45?
this is a part of the code I am using with the attacched workspace that has to be loaded by the code to run:
% this is one of many funcions in my code
% each load step can be either dynamic or static over a fictitius time
% interval. For now I have a linearised and non-linear g-alpha solvers validated against abaqus,
% but I wanted to use also a solver with a doubled dof system and use ode45
% though I always struggle to get the right formatting to use embedded
% This is an extrapolated demo of my code and for now it won't
% work for multple dynamic load steps solved with ode45.
% This demo is thought to help the kind person that is willing to help me and therefore it will stop at the end of te first load steps that should be required
% the problematic part ishighlighted below...
clc
clear all
close all
t0=0;
Step{1}=RefConfig;
d0 = zeros(length(UnknownDispIndex),1); % initial displacement is zero
v0 = d0*0;
a0 = d0*0;
dMtot=[];
vMtot=[];
aMtot=[];
step_counter=1;% used in the function to moove among different TIME steps.. independently from what type of step if linear or dynamic
timeline=0;
for i=1:input.load_step_number-1 % there are 2 steps in the analysis but I am just interested in getting the first one running
% my problem is with even one step of ode45 starting from line 50
if input.StepSolverType{i}==1 % linear static solver
%other solver not in this demo

elseif input.StepSolverType{i}==2 % for now only ninear but maybe one day also non linear but on th full model
% g_alpha solver
if input.AlgorithmType{1}==1
if input.reduced{i}==1
err('Galpha can only be used for full models for now')
else
if input.discreteVScontinuous==1
%other solver not in this demo
else % discreteVScontinuous==2 ... for continuously varing force vector with time

if i==1 % just the very first time step require this update of glob force vec the others are fine as they are updated in the loop
else
end
%other solver not in this demo
end
end
% ABC matrix with ode45

elseif input.AlgorithmType{2}==2 % ABC matrix with ode45
if input.reduced{i}==1
err('I still have to do all the projections on the eigenvector matrix')
else
n=length(d0);
A = [zeros(n,n), eye(n);-Step{step_counter}.M22\Step{step_counter}.K22, zeros(n,n)];
if input.discreteVScontinuous==1 % discrete force matrices
err('I do not think it is worth to spend time on a staet vector augmented solver for forces given as fixxed arays in time and do not take advantage of ode45 otherwise just use galpha and do not ask too much')
else % discreteVScontinuous==2 ... for continuously varing force vector with time
state_vec0=[d0;v0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the problematic part %%%%%%%%%% I can write a
% generalized alpha solver but I am struggling
% with running embedded matlab functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=@(state_vec,time)eqfunODE(n,A,Step{step_counter}.M22,Step{step_counter}.K22,state_vec,time,TotalGlobForceVecReduced{i});
[t,Y001] = ode45(@eqfunODE,[input.time_check_point{i},input.time_check_point{i+1}],state_vec0);
[t,Y002] = ode45(f,[input.time_check_point{i},input.time_check_point{i+1}],state_vec0); %differnet way of writing just to check if it is the same and it is
end
end
else
end
elseif input.StepSolverType{i}==3 % I might put a non linear static solver with geometrical stiff
else
end
% use for continuity over loadsteps
% t0=tc;
% d0 = dM(:,end);
% v0 = vM(:,end);
% a0 = aM(:,end);
% dMtot=cat(2,dMtot,dM);
% vMtot=cat(2,vMtot,vM);
% aMtot=cat(2,aMtot,aM);
end
% this uually are the output of this function
% sol.dMtot=dMtot;
% sol.vMtot=vMtot;
% sol.aMtot=aMtot;
% sol.Step=Step;
% sol.timeline=timeline;
function [dydt]=eqfunODE(n,A,M22,K22,state_vec,time,TotalGlobForceVecReduced)
y1=zeros(n,1);
y2=zeros(n,1);
dydt_1=zeros(n,1);
dydt_2=zeros(n,1);
dydt=zeros(2*n,1);
y1=state_vec(1:n);
y2=state_vec(n+1:2*n);
dydt_1=y2;
dydt_2=-M22\K22*y1+M22\TotalGlobForceVecReduced(time);
dydt = [dydt_1; dydt_2];
% other approach not working
% Fexternal=[zeros(n,1);M22\TotalGlobForceVecReduced(time)];
%

%
%
% equdif=A*state_vec(1:2*n,1)+Fexternal(time);
end