how to solve 2nd order coupled system of differential equations with heaviside function using ode45 solver?

i was trying to solve a system of coupled differential equations using ode45 solver with one of the equations having a heaviside function but i guess that it is assuming the value of zero of the heaviside function. please help.

syms P k1 k0 mw Izw b k22 mt Iyg kt dg ct Ixg v w kf delta;
P = 253;
k1 = 1.8;
k0 = 0.025;
mw = 2.53;
Izw = 0.0055;
b = 0.049;
k22 = 171;
mt = 0.41;
Iyg = 0.00013;
kt = 0.0045;
kf = 0.45;
delta = 0.009; 
dg = 1.45/2;
ct = 0.00063;
Ixg = 0.00037;
v = 100;
w = 10.47*30;
    couplode = @(t,y) [y(2); -((2*k22/v*(mt+mw))*y(2))+((2*k22/(mt+mw))*y(3))-((2*P*k1/(mt+mw))*y(1))-(kf*(y(1)+delta)*heaviside(y(1)-delta));y(4);(-2*(P*k0*b/(Izw+Iyg))*y(3))+((Ixg*w/(Izw+Iyg))*y(6));y(6);((-kt*dg^2/Iyg)*y(5))-((Ixg*w/Iyg)*y(4))-((ct*dg^2/Iyg)*y(6))];
    [t,y] = ode45(couplode, [0 1000*pi], [0.001;0;0.314;0;0;0]); 
    figure(1)
    plot(t, y)
    hold on;
grid;

.pl

Answer 1

[Star Strider]

Numerical ODE solvers do not do well across non-differentiable discontinuities. The heaviside function in MATLAB is differentiable, so you likely can use it if you want.

I do not see where you have called heaviside or anything similar to it. What do you want to do?

[Ronak Prateek]

couplode = @(t,y) [y(2); -((2*k22/v*(mt+mw))*y(2))+((2*k22/(mt+mw))*y(3))-((2*P*k1/(mt+mw))*y(1))-(kf*(y(1)+delta)*heaviside(y(1)-delta));y(4);(-2*(P*k0*b/(Izw+Iyg))*y(3))+((Ixg*w/(Izw+Iyg))*y(6));y(6);((-kt*dg^2/Iyg)*y(5))-((Ixg*w/Iyg)*y(4))-((ct*dg^2/Iyg)*y(6))];

[Ronak Prateek]

i wanted to obtain time history of the system. but i believe the heaviside function takes 0 value.

[Star Strider]

I did not see the heaviside call before.

The system runs without error and appears to produce appropriate output.

If you want to see what the heaviside function is actually doing, add a specific plot call for it:

    figure(1)
    plot(t, y)
    hold on
    plot(t, heaviside(y(:,1)-delta), '--k')
    hold off

Here, it plots the result as a dashed black line.

[Ronak Prateek]

this is the system of equation which i was trying to solve with heaviside functions being highlighted below. my code displays the same result with or without the consideration of heaviside function terms. i believe the first highlighted term as shown below can be represented in the matlab by

(kf*(y(1)+delta)*heaviside(-y(1)-delta))

is it correct or not? please tell. sorry for not being clear.