## Problem in simulation ODE45 Index exceeds array bounds.

### Fernando Osorio (view profile)

on 10 Nov 2018
Latest activity Answered by Fernando Osorio

on 11 Nov 2018

### Fernando Osorio (view profile)

Hi. I am new in Matlab and I am trying to solve the following differential equation.
function [xdot] = HW5P3(t,x)% x is theta
%Function equations
v=12.*sin(2.*pi.*(1./2).*t); %sine source
L=1-(2./pi).*abs(x(1));
hdot=-((1./pi).*(sign(x(1))./(L.^2)).*(x(4)).^2)-10.*x(3)-x(1);
end
[t,x]=ode45('HW5P3',[0 5],[1,0,0]);

R2018a

### Fernando Osorio (view profile)

on 11 Nov 2018

I added more elements to the code and also new conditions but this one worked for me.
% Function
function [xdot y] = HW5P3(t,x)% x is theta
%Function equations
global a b B J K vo R Lo thetaout
theta=x(1);
h=x(2);
lambda=x(3);
v=(t<0.5).*(0)+((t>=0.5)&(t<4.5)).*(vo*sin(2*pi*(1/0.1)*t))+(t>=4.5).*(0); %Voltage source starts at t=0.5sec
L=(theta<-thetaout).*(Lo)+((theta>=-thetaout)&(theta<thetaout)).*(a+b.*abs(theta))+(theta>=thetaout).*(Lo); %Inductance change with theta
Ldot=b*sign(theta);%Derivative of L
hdot=((Ldot/(2*(L)^2))*lambda^2)-((B/J)*h)-K*theta;
y=[L;v;Ldot];
end
% code
global a b B J K vo R Lo thetaout
a=1;
b=-2/pi;
B=0.01;%Damping factor for bearings
J=0.001;%Inertia of the shaft
K=1;%Rotational Spring of the shaft
Lo=0.1;%minimum inductance when the iron of the blade is outside the gap
vo=120;%source voltage pk-pk
R=10; %The resistance value was to high and the source wasn't giving enough current. I changed it to 10ohms
thetaout=1.4137;%limit of the blade to work by the effect of L(theta) out of this angle L=Lo
lambdai=0;%Initial condition
thetai=1.4137;%Initial COndition
hi=0;%Initial COndition
[t,x]=ode45('HW5P3',[0 5],[thetai,hi,lambdai]); %calling of function
L=length(t);
N = length(t);
for i=1:N %getting values of output from the function for every time point
[xd,y] = HW5P3(t(i),x(i,:));
Ln(i) = y(1);%L(theta)
vs(i) = y(2);%Voltage source
Lndot(i) = y(3);%Derivative of L(theta)
end
theta=x(:,1);
h=x(:,2);
lambda=x(:,3);
is=lambda/L;
Trotspring=K*theta;
figure(1);
subplot(3,1,1)
subplot(3,1,2)