function xdot=My2ndAugModel(t,x) x1=x(1); % Angle
x2=x(2); % Velocity
x3=x(3); % Ft
x4=x(4); % Lm
x5=x(5); % Lmdot
m=75;
g=9.8;
k=0.25;
J=2.5;
Mm=0.435;
Fmc=69.45;
lmc=0.125; Kt=63000; rp=0.05; Km=6700; Bm=257.1;
k1=1.5;
k2=0.9;
k3=2.5;
k4=0.5;
Bp=5.57;
thetamin=(pi/180);
thetamax=5*(pi/180);
Mp=(k1*exp(-k2*(x1-thetamin))-k3.*exp(-k4*(thetamax-x1))-Bp*x2);
x1dot=x2;
x2dot=(((m*g*k)*cos(x1))-(rp*x3)+Mp)/J;
x3dot=Kt*(rp*x2-x5);
x4dot=x5;
x5dot=(x3-((Km*(x4-lmc)+Fmc)+Bm*x5))/Mm;
xdot=[x1dot;x2dot;x3dot;x4dot;x5dot];
end %%%%%%%%%%%%%%%% With Run File %%%%%%%%%%%%%%%%%%%% clc
close all
clear all
t=linspace(0,3,500);
F0=1200;
xinit=[0 0 0 0.0959 0];
options=odeset('RelTol',1e-4,'AbsTol',1e-5);
%[t,x]=ode45(@My2ndAugModel,t,xinit,options);
[t,x]=ode45(@(t,x) My2ndAugModel(t,x),t,xinit,options)
angle=x(:,1)*(180/pi);
vel=x(:,2);
Ft=x(:,3);
Lm=x(:,4);
Lmdot=x(:,5);
subplot(321);
plot(t,angle);ylabel('Angle (deg)')
subplot (322)
plot(t,vel);ylabel('\thetadot (rad/sec)')
subplot(323)
plot(t,Ft);ylabel('Ft')
subplot(324)
plot(t,Lm);ylabel('Lm')
subplot(325)
plot(t,Lmdot);ylabel('Lmdot')
