A function for getting F and My from simulation results

FBmodel.m

function dydt=FBmodel(t,y)
%global y2

% Model file to be initiated by FBstart.m
% extract variable from y-vector

X=y(1);
S=y(2);
V=y(3);

% Constants

qSmax=1.6;
Ks=0.1;
qm=0.04;
Yem=0.5;
Si=500;
F0=0.03;
SFR=0.3;
Fmax=0.8;

% algorithm

F=F0*exp(SFR*t);
if F>Fmax
    F=Fmax;
end

qS=qSmax*S/(S+Ks);
My=(qS-qm)*Yem;
dXdt=-F/V*X+My*X;
dSdt=F/V*(Si-S)-qS*X;
dVdt=F;

% make a dydt column vector

dydt=[dXdt; dSdt; dVdt];

% Store non-differential variables in y2

%y2=[y2;[t,F,My]];

fbstart.m

% FBstart; Initiation file for fed-batch simulation.
% Requires a seperate model file

clear all

global y2
y2=[];                  % for storage of non differential equation variables
tspan=[0 50];           % % time scale  

% enter initial values and locate in column vector

X=0.5;
S=0.1;
V=40;

y=[X;S;V;];

% call ODE solver and the model file

[t y]=ODE23s('FBmodel',tspan,y);

[simulated_My, simulated_F] = GetValues(t,y); % Get the simulated My and F values from this run

[t_prev, y_prev] = textread('') %do whatever you need here to load in the T and Y values from the text-file
[prev_My prev_F] = GetValues(t_prev, y_prev); %get My and F values from the text-data

prev_My = interp1(t_prev, prev_my, t);
prev_F = interp1(t_prev, prev_F, t);

% This stuff is no longer needed
%if isempty(y2)==0    
%    % eliminates duplicates
%    
%    y2(find(diff(y2(:,1))<diff(tspan)/1000),:)=[];
%    y2(find(diff(y2(:,1))<0),:)=[];
%    
%    % match to y-vector size
%   
%    y2=interp1(y2(:,1),y2(:,2:length(y2(1,:))),t);
%   
%    % merge with y-vector
%    
%    y=[y,y2];
% end

% scale max values for scaling in graph

ymax=[100,1,100,1,1];   %X,S,V,F, My

% scale values to a 0-100 scale

for i=1:length(ymax)
    yscaled(:,i)=y(:,i)/ymax(i)*100;
end

% plot and label

yplot=plot(t,yscaled);
set(gca,'YLim',[0 100])
legend('X: 0-100 g/L','S: 0-1 g/L','V: 0-100 L','F:0-1 L/h','My: 0-1 /h')
xlabel('time (hrs)')
title('Fed-batch with eponential/constant feed')
figure(gcf)

Getvalues.m

function [F, My] = GetValues(in_t, indata)

X=indata(:,1);
S=indata(:,2);
V=indata(:,3);

% Constants

qSmax=1.6;
Ks=0.1;
qm=0.04;
Yem=0.5;
Si=500;
F0=0.03;
SFR=0.3;
Fmax=0.8;

% algorithm

F=min(F0*exp(SFR*in_t), Fmax);
if F>Fmax
    F=Fmax;
end

qS=qSmax*S/(S+Ks);
My=(qS-qm)*Yem;

% make a dydt column vector


% Store non-differential variables in y2