FuelCellVoltAmpSolver.m

clear;
clc;

Temp = 477;
power_load = 2198;

A = 0.023951368224792 * Temp + 23.9241562583015;
B = 0.003480859912024 * Temp - 2.19986938582928;
C = -0.0001779207513 * Temp + 0.104916556604259;
D = 5.0656524872309E-06 * Temp - 0.002885372247954;
E = -6.42229870072935E-08 * Temp + 3.58599071612147E-05;
F = 3.02098031429142E-10 * Temp - 1.66275376548748E-07;

Volts = 23.583; % <------- ADJUST TO THE LAST VOLTAGE PRINTED ON THE LAST RUN
Amperes = 66.177; % <---- ADJUST TO THE LAST CURRENT PRINTED ON THE LAST RUN

Vlast = 0;
Ilast = 0;
loadResistance = (784) / (power_load);

v_data = [];
i_data = [];
step_data = [];
ii = 1;
    
while true
  v_data(ii) = Volts;
  i_data(ii) = Amperes;
  step_data(ii) = ii;
  
  Vlast = Volts;
  Ilast = Amperes; 	
  
  Amperes = Volts / (loadResistance)
  Volts = A +...
  B * Amperes +...
  C * Amperes*Amperes +...
  D * Amperes*Amperes*Amperes +...
  E * Amperes*Amperes*Amperes*Amperes +...
  F * Amperes*Amperes*Amperes*Amperes*Amperes
  
  
  if(abs(Volts-Vlast)<0.001)
    ii++;
    v_data(ii) = Volts;
    i_data(ii) = Amperes;
    step_data(ii) = ii;
    break;
  endif
  
  ii++;
  if(ii>100)
    break;
  endif
endwhile

[ax, plt1, plt2] = plotyy(step_data,v_data,step_data,i_data);
grid on;
grid minor;
ylabel(ax(1),"Voltage");
ylabel(ax(2),"Current");
xlabel("Iteration")

P = Volts*Amperes