clear all; close all;


R=20; L=2e-3; C=40e-6;
global E=40;
A=[-1/(R*C) 1/C; -1/L 0];
b1=[0 0]'; b2=[0 E/L]';
global u_ref=32;
x_1=[u_ref u_ref/R]';

global  t_off=2.52683923950417e-04;
global  t_on=3.05171993781422e-04;


%Equacion diferencial
function retval=f2(x,u,t,params)
  retval=params.A*x + (params.b2).*u(t);
end

function retval=control_opt(t)
  freq=1e7;
  global t_off;
  global t_on;
  global E;
  global u_ref;
  
  
  if t<=t_on
    retval=1;
    return;
  end
  if t<=t_on+t_off
    retval=0;
    return;
  end
  retval=(square((t - t_on-t_off)*freq/(2*pi), 0.8)+1)/2; % 2pi para la normalizacion, see help square  
end

params=struct("A", A, "b1", b1, "b2", b2);

t=[linspace(0,t_on-1/100*t_on,100) linspace(t_on,t_on+t_off-1/100*t_off,100) linspace(t_on+t_off,1e-3+t_on+t_off,10000+1)];

%The bug is triggered here. To see the difference uncomment the line with ode4
%and compare the results.
x=lsode(@(x,tau) f2(x, @control_opt ,tau, params), [0 0]', t);

%x=ode4(@(tau,x) f2(x, @control_opt ,tau, params), t, [0 0]');


plot(t, x(:,1),t,x(:,2), t, arrayfun(@control_opt ,t))
