function xdot=plant(t,x);
% PLANT -   function to be used with ode23s
%
% INPUT
%   t   -   time vector according to the time span specified
%   x   -   values of x to be simulated
%
% OUTPUT
%   xdot    -   time derivative of the plant
global R11 R12 R21 R22 C1 C2 r1 r2 I1 I2 theta
xdot =[((x(1)/R11) + (x(2)/R12) + I1 - (sigInv(x(1), theta)/r1))./(C1.*dsigInvdx(x(1), theta));
       ((x(1)/R21) + (x(2)/R22) + I2 - (sigInv(x(2), theta)/r2))./(C2.*dsigInvdx(x(1), theta))];
   
function u = sigInv(x, theta)
u = theta*(log(x) - log(1-x));

function y = dsigInvdx(x, theta)
y = (theta./x) + theta./(1-x);