function [Psires,dPsires_dxdot,dPsires_dx,dPsires_du,dPsires_dd,dPsires_dw] ...
    = psi_equation(x,xdot,u,geop,~,trap,BCs,model)

nx = model.dims.nx;
nu = model.dims.nu;
nrg = model.rgrid.nrhogauss;

%%% Physics
% General equation structure
% 0 = - mz*dz/dt + d/drho*dz*dz/drho + az*dz/drho + hz*z + lz + sz
% f = - Mz*zhatdot - Dz*zhat + Az*zhat + Hz*zhat + lz + sz

% Physical constants
mu0 = 4e-7*pi;

% combine all psi sources
jniB = trap.jbsB + trap.jauxB + trap.jreB + trap.jdistB; % <jni dot B>

BCval = BCs.psiValue;

mpsi_varying  = geop.Phib.^2*(trap.signeo./geop.F.^2);
apsi_varying  = geop.Phibdot.*geop.Phib.*trap.signeo./geop.F.^2;
dpsi_varying  = geop.g23or;
hpsi_varying  = zeros(nrg,1);
spsi_varying  = geop.Phib*geop.Vp./geop.F.^2.*jniB;
bcpsi_varying = geop.Phib/geop.Fb * BCval;
     

%%% Matrices
% Mpsi (function of Te)
Mpsi = zeros(model.psi.nsp);
Mpsi(model.mats.Mpsi.XIindex) = model.mats.Mpsi.XI*mpsi_varying;

% Apsi
Apsi = zeros(model.psi.nsp);
Apsi(model.mats.Apsi.XIindex) = model.mats.Apsi.XI*apsi_varying;

% Dpsi
Dpsi = zeros(model.psi.nsp,model.psi.nsp);
Dpsi(model.mats.Dpsi.XIindex) = model.mats.Dpsi.XI*dpsi_varying;

% Hpsi
Hpsi = zeros(model.psi.nsp);
Hpsi(model.mats.Hpsi.XIindex) = model.mats.Hpsi.XI*hpsi_varying;

KKpsi = Apsi + Hpsi - Dpsi; % total term for contributions proportional to Psihat

% boundary condition
bcpsi_constant = zeros(model.psi.nsp,1);
bcpsi_constant(end) = 16*pi^3*mu0/model.psi.scal;
lPsi = bcpsi_constant*bcpsi_varying; % influence on inductive source through b.c. (eq. 19 in equation summary)

% sources
sPsi = model.mats.Spsi.XI*spsi_varying;

% evaluate Psi equation residual ("F" in FF thesis eq7.40)
Psihatdot = xdot(model.psi.xind);
Psihat    = x(model.psi.xind);

%%
Psires = -Mpsi*Psihatdot + KKpsi*Psihat + lPsi + sPsi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Jacobian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargout > 2;
    % general structure:
    % f = - M*psidot + A*psi - D*psi + H*psi + s + BC
    % df/dxdot = -dM/dxdot * psidot - M*(dpsidot/dx) + dA/dxdot * psi           - dD/dxdot * psi          + dH/dxdot *psi            + ds/dxdot
    % df/dx    = -dM/dx    * psidot                  + dA/dx * psi + A(dpsi/dx) - dD/dx * psi+ D(dpsi/dx) + dH/dx * psi + H(dpsi/dx) + ds/dx
    % df/du    = -dM/du    * psidot                  + dA/du * psi              - dD/du * psi             + dH/du * psi              + ds/du 
    
    dPsi_dx = zeros(model.psi.nsp,nx);
    dPsi_dx(1:model.psi.nsp,model.psi.xind) = eye(model.psi.nsp);
    
    %% Psi equation
    
    % Sources
    djniB_dx    = trap.djauxB_dx + trap.djbsB_dx;
    djniB_du    = trap.djauxB_du; % + 0    ;
    djdisB_dd   = [model.dist.psi.Lamgauss*model.equi.B0,zeros(model.rgrid.nrhogauss,numel(model.dist.te.ind))]; 
    djdisB_dw   = [model.dist.psi.Lamgaussw*model.equi.B0,zeros(size(model.dist.psi.Lamgaussw))]; 
    
    iFg2 = 1./geop.F.^2;
    PhibVpiFg2 = geop.Phib*(geop.Vp.*iFg2);
    dmpsivarying_dx = bsxfun(@times,geop.Phib.^2.*iFg2,trap.dsigneo_dx);
    dapsivarying_dx = bsxfun(@times,geop.Phibdot.*geop.Phib.*iFg2,trap.dsigneo_dx);
    dspsivarying_dx = bsxfun(@times,PhibVpiFg2,djniB_dx);
    dspsivarying_du = bsxfun(@times,PhibVpiFg2,djniB_du);
    dspsivarying_dd = bsxfun(@times,PhibVpiFg2,djdisB_dd);
    dspsivarying_dw = bsxfun(@times,PhibVpiFg2,djdisB_dw);
    dlPsi_dx  = 0; %BC does not depend on x
    dlPsi_du = geop.Phib/geop.Fb*(bcpsi_constant*BCs.dpsiValue_du); % BC depends on u (see above)
    
    
    % M matrices
    zx = zeros(model.psi.nsp,nx);
    zu = zeros(model.psi.nsp,numel(u));
    dMpsi_dxdot_Psihatdot = zx; % M does not depend on xdot
    
    dMpsi_dx_Psihatdot = zeros(model.psi.nsp,nx); % init
    dMpsi_dxi = zeros(model.psi.nsp); % init
    
    if any(any(dmpsivarying_dx))
        dMpsi_dxi_entries = model.mats.Mpsi.XI*dmpsivarying_dx;
        for ii=1:nx;
            dMpsi_dxi(model.mats.Mpsi.XIindex) = dMpsi_dxi_entries(:,ii);
            dMpsi_dx_Psihatdot(:,ii) = dMpsi_dxi*Psihatdot;
        end
    end
    
    %%
    dMpsi_du_Psihatdot = zu; % Mpsi does not depend on u
    
    % A matrices
    dApsi_dxdot_Psihat = zx; 
    
    dApsi_dx_Psihat = zeros(model.psi.nsp,nx); % init
    dApsi_dxi = zeros(model.psi.nsp); % init
    if any(any(dapsivarying_dx))
        dApsi_dxi_entries = model.mats.Apsi.XI*dapsivarying_dx;
        for ii=1:nx;
            dApsi_dxi(model.mats.Apsi.XIindex) = dApsi_dxi_entries(:,ii);
            dApsi_dx_Psihat(:,ii) = dApsi_dxi*Psihat;
        end
    end
    dApsi_du_Psihat = zu;
    
    % D matrices
    dDpsi_dxdot_Psihat = zx; % Dpsi does not depend on xdot
    dDpsi_dx_Psihat    = zx; % Dpsi does not depend on x
    dDpsi_du_Psihat    = zu; % Dpsi does not depend on u
    
    % H matrices
    dHpsi_dxdot_Psihat = zx; 
    dHpsi_dx_Psihat    = zx;
    dHpsi_du_Psihat    = zu;
    
    % Spsi
    dSpsi_dxdot = zx; % Spsi does not depend on xdot
    dSpsi_dx    = model.mats.Spsi.XI*dspsivarying_dx + dlPsi_dx;
    dSpsi_du    = model.mats.Spsi.XI*dspsivarying_du + dlPsi_du;
    dSpsi_dd    = model.mats.Spsi.XI*dspsivarying_dd;
    dSpsi_dw    = model.mats.Spsi.XI*dspsivarying_dw;

    % Jacobians for Psi equation
    dPsires_dxdot = -dMpsi_dxdot_Psihatdot - Mpsi*dPsi_dx + dApsi_dxdot_Psihat- dDpsi_dxdot_Psihat + dHpsi_dxdot_Psihat + dSpsi_dxdot;
    dPsires_dx    = -dMpsi_dx_Psihatdot    + dApsi_dx_Psihat + Apsi*dPsi_dx - Dpsi*dPsi_dx - dDpsi_dx_Psihat + dHpsi_dx_Psihat + Hpsi*dPsi_dx   + dSpsi_dx;
    dPsires_du    = -dMpsi_du_Psihatdot       + dApsi_du_Psihat            - dDpsi_du_Psihat   + dHpsi_du_Psihat + dSpsi_du;
    dPsires_dd    = dSpsi_dd;
    dPsires_dw    = dSpsi_dw;
      
end

return