% 2 Degree of Freedom (With Damping)
% Dr Owen Williams
% Institute of Science and Art
% University of Wales Trinity Saint David

clear; clc
syms x(t) phi(t) Y
Dx = diff(x) 
D2x = diff(x,2) 
Dphi = diff(phi) 
D2phi = diff(phi,2)
m=200; kf=12000; kr=12000; cf=1400; cr=1400; I=28; p=1.4; b=0.65 


eqn1 = m*D2x+(cf+cr)*Dx(t)+(cf*(p-b)-cr*b)*Dphi(t)+(kf+kr)*x(t)+(kf*(p-b)-kr*b)*phi(t)==0; 
eqn2 = I*D2phi+(cf*(p-b)-cr*b)*Dx(t)+(cf*(p-b)^2+cr*b^2)*Dphi(t)+(kf*(p-b)-kr*b)*x(t)+(kf*(p-b)^2+kr*b^2)*phi(t)==0;

[V] = odeToVectorField(eqn1,eqn2)
vars = x(t)
M = matlabFunction(V,'vars',{'t', 'Y'})
interval = [0 2]
xSol=ode45(M,interval,[1 0 0 0])

tValues = linspace(interval(1),interval(2),1000);
yValues = deval(xSol,tValues,1); %number 1 denotes first solution likewise you can mention 2 for the next solution
plot(tValues,yValues)
hold on
yValues = deval(xSol,tValues,3); 
plot(tValues,yValues,'r')