// Periodic input to ODEs
//
// 1 from -Pi to 0 and 0 for 0 to Pi
// Square Wave
function fout = sqwave(t)
	x = pmodulo(t,2*%pi);
	fout = zeros(x);
	fout(x>%pi) = 1;
endfunction
// Triangle down first
function fout = tdownwave(t)
	fout = abs(pmodulo(t,2*%pi)-%pi);
endfunction
// Triangle up first
function fout = tupwave(t)
	fout = %pi-tdownwave(t);
endfunction
// Saw Tooth
function fout = sawwave(t)
	fout = (%pi-pmodulo(t,2*%pi))/ 2;
endfunction

scf(0);
clf()
t=0:0.05:10*%pi;
plot(t,sqwave(t),'r-')
plot(t,tupwave(t),'g-')
plot(t,tdownwave(t),'b-')
plot(t,sawwave(t),'k-')

function x = force(t)
	x = tdownwave(t)
endfunction

// The ODE as a system y'' + 4y = force(t) is
// w = [y, z] where y' = z and z' = y'' = -4 y + force(t)
function wdot = f(t, w)
wdot = [w(2,1); -4*w(1,1) + force(t); ]
endfunction

w=ode([0;0],0,t,f);
scf(1);
clf();
plot(t,force(t),"b",t,w(1,:),"r");
xtitle("UnDamped Input in Blue, Output in Red","time t","y");


// The Damped ODE system 
// u'' + 0.125 u' + 4 u = force(t)
// u' = z
function wdot = f(t,w)
wdot = [w(2,1); -0.125*w(2,1) - 4*w(1,1) + force(t); ];
endfunction

w=ode([0;0],0,t,f);
scf(2);
clf();
plot(t,force(t),"b",t,w(1,:),"r");
xtitle("Damped Input in Blue, Output in Red","time t","y");