paso=0.5;
tiempo=0:paso:250;
function dhdt=dhdt(h,t)
    if t<=30
        f=3;
    else
        f=3+0.01*t;
        if f>=5;
              f=5;
        end
    end
    global x
    dhdt=f-x*h;
endfunction %%las funciones unicamente pueden ver las varibales que están definidas adentro de ellas. En la función se toma un valor de x, sin embargo se define fuera de la funcion y no la puede reconocer. Por lo tanto, se tiene que poner como variable global. Se tiene que definir adentro y afuera de la función.
%%la ecuación del controlador derivativo/integral es la siguiente: pv=bias+kc(er+ti*int(Edt)+td*dE/dt)
ho=1.5;
h=h0;
sp=6; %%set point: el valor que alcanza el estado estable sin perturbación%%
kc=-0.3;
ti=0;
td=0;
inter=0;
global x
for i=1:length(tiempo)-1 %%se mueve a lo largo del vector y a lo largo de este, se establecen los puntos en el tiempo en los que se resuelven las ecuaciones diferenciales hasta completar el intervalo a analizar. En cada intervalo de tiempo se puede calcular como varía el error con respecto al tiempo%%
    if tiempo(i)<=30;
    er(i)=0;
    dedt(i)=0;
    x=0.5;   %%durante los primeros 30 segundos se establece un valor del error de 0 porque en ese intérvalo no hay perturbaciones
    else;
        er(i)=sp-h(i);
        dedt(i)=(er(i)-er(i-1))/paso; %%se calcula el error con respecto al tiempo en un intervalo dado
        x=0.5+kc*(er(i-1)+ti*inter+td*dedt(i));%%en este paso empiezan a funcionar los controladores
    end
    inter=inter+er(i)*paso;
    h(i+1)=lsode("dhdt",ho,[tiempo(i),tiempo(i+1)])(2); %%se generarían los siguientes valores de "h"
    %%la ecuación diferencial da dos valores: el valor inicial y el valor final...en este comando se le pide al programa que tome el   segundo valor del vector
    h0=h(i+1);
end
plot(tiempo,h)