function simple_solver

%% solve system of two ODEs:
% dm/dt = f(p) - d1*m
% dp/dt = k2*m - d2*p
% f(p) is either a decreasing or increasing Hill function

    %% set model parameters
    
    d1 = 0.1; %m degradation rate
    d2 = 0.08; %p degradation rate
    k1 = 1; %max transcription rate
    k2 = 1; %translation rate
    p0 = 50; %Hill function threshold
    n = 1; %Hill function steepness
    a = 1; %a=0 is positive feedback; a=1 is negative feedback
           % a must be set to either 0 or 1 (not any other value)

    if (a~=0)&(a~=1)==1 
        error('incorrect value of a; must be 0 or 1')
    end
    
    %% set time interval and initial values
    
    tstart = 0;
    tend = 300;
    
    m_0 = 1;
    p_0 = 0;

    %% solve the ODEs -- do not change
    
    sol=ode45(@odes,[tstart,tend],[m_0;p_0]);
    tout=tstart:1:tend;
    stout=deval(sol,tout);

    %% plot the solution
    
    figure(1), clf
    subplot(2,1,1)
    set(gca,'FontSize',16)
    plot(tout,stout,'LineWidth',2)
    xlabel('t')
    ylabel('amount of m and p')
    legend('m(t)','p(t)')

    
    %% plot the nullclines
    
    subplot(2,1,2)
    set(gca,'FontSize',16)
    x=linspace(0,4*p0,100);
    plot(x,f(x)/d1,'b',x,d2*x/k2,'r',stout(2,:),stout(1,:),'k','LineWidth',2)
    xlabel('p'), ylabel('m')
    if a==0
    legend('m nullcline', 'p nullcline','trajectory','Location','NorthWest')
    else
    legend('m nullcline', 'p nullcline','trajectory','Location','NorthEast')
    end
  
    
    %% function to specify the ODEs
    
    function dydt=odes(t,y)
        m=y(1);
        p=y(2);
        dydt=[f(p)-d1*m; k2*m-d2*p];
    end

    %% function to evaluate function f(p)
    
    function out=f(p)
        out=a*k1./(1+(p./p0).^n)+(1-a)*k1*(p.^n)./(p0.^n + p.^n);
    end
end