Unrecognized function or variable 'RK'.

% Example5_5.m
% Solution to Example 5.5. This program calculates and plots
% the concentration of cell mass, concentration of penicillin,
% optimal temperature profile, and adjoint variable of a batch
% penicillin fermentor.  It uses the function COLLOCATION to 
% solve the set of system and adjoint equations.
clear
clc
clf
% Input data
w = input(' Enter w''s as a vector : ');
y0 = input(' Vector of known initial conditions = ');
yf = input(' Vector of final conditions = ');
guess = input(' Vector of guessed initial conditions = ');
fname = input('\n M-file containing the set of differential equations : ');
fth = input(' M-file containing the necessary condition function : ');
n = input(' Number of internal collocation points = ');
rho = input(' Relaxation factor = ');
% Solution of the set of differential equations
[t,y] = collocation(fname,0,1,y0,yf,guess,n,rho,[],w,fth);
% Temperature changes
for k = 1:n+2
    options=optimset;
    theta(k) = fzero(fth,30,options,y(:,k),w);
end
% Plotting the results
subplot(2,2,1), plot(t,y(1,:))
xlabel('Time')
ylabel('Cell')
title('(a)')
subplot(2,2,2), plot(t,y(2,:))
xlabel('Time')
ylabel('Penicillin')
title('(b)')
subplot(2,2,3), plot(t,y(3,:))
xlabel('Time')
ylabel('First Adjoint')
title('(c)')
subplot(2,2,4), plot(t,theta)
xlabel('Time')
ylabel('Temperature (deg C)')
title('(d)')
function f = Ex5_5_func(t,y,w,fth)
% Function Ex5_5_func.M
% This function introduces the set of ordinary differential
% equations used in Example 5.5.
% Temperature
options=optimset;
theta = fzero(fth,30,options,y,w);
% Calculating the b's
b1 = w(1) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
if b1<0, b1=0; end
b2 = w(4) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
if b2<0, b2=1e-6; end
b3 = w(5) * (1-w(2)*(theta-w(6))^2) / (1-w(2)*(25-w(6))^2);
if b3<0, b3=0; end
% Evaluating the function values
f(1) = b1*y(1) - b1/b2*y(1)^2;
f(2) = b3*y(1);
f(3) = -b1*y(3) + 2*b1/b2*y(1)*y(3) - b3;
f = f';	% Make it a column vector
function ftheta = Ex5_5_theta(theta,y,w)
% Function Ex5_5_theta.M
% This function calculates the value of the necessary condition
% as a function of the temperature (theta). It is used in solving
% Example 5.5.
% Calculating the b's
b1 = w(1) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
db1 = w(1)*(-w(2))*2*(theta-w(3)) / (1-w(2)*(25-w(3))^2);
b2 = w(4) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
db2 = w(4)*(-w(2))*2*(theta-w(3)) / (1-w(2)*(25-w(3))^2);
b3 = w(5) * (1-w(2)*(theta-w(6))^2) / (1-w(2)*(25-w(6))^2);
db3 = w(5)*(-w(2))*2*(theta-w(6)) / (1-w(2)*(25-w(6))^2);
% The function
ftheta = y(3)*(y(1)*db1-y(1)^2*(db1*b2-db2*b1)/b2^2)+y(1)*db3;
function [x,y] = collocation(ODEfile,x0,xf,y0,yf,guess,n,rho,tol,varargin)
%COLLOCATION Solves a boundary value set of ordinary differential
% equations by the orthogonal collocation method.
%
%   [X,Y]=COLLOCATION('F',X0,XF,Y0,YF,GAMMA,N) integrates the set of
%   ordinary differential equations from X0 to XF by the Nth-degree
%   orthogonal collocation method.  The equations are contained in
%   the M-file F.M.  Y0, YF, and GAMMA are the vectors of initial
%   conditions, final conditions, and starting guesses respectively.
%   The function returns the independent variable in the vector X
%   and the set of dependent variables in the matrix Y.
%	
%   [X,Y]=COLLOCATION('F',X0,XF,Y0,YF,GAMMA,N,RHO,TOL,P1,P2,...)
%   uses relaxation factor RHO and tolerance TOL for convergence
%   test.  Additional parameters P1, P2, ... are passed directly to
%   the function F.  Pass an empty matrix for RHO or TOL to use the
%   default value.
%
%   See also SHOOTING
% (c) N. Mostoufi & A. Constantinides
% January 1, 1999
% Initialization
if nargin < 7 | isempty(n)
    n = 1;
end
if nargin < 8 | isempty(rho)
    rho = 1;
end
if nargin < 9 | isempty(tol)
    tol = 1e-6;
end
y0 = (y0(:).')';			% Make sure it's a column vector
yf = (yf(:).')';			% Make sure it's a column vector
guess = (guess(:).')';	% Make sure it's a column vector
% Checking the number of guesses
if length(yf) ~= length(guess)
    error(' The number of guessed conditions is not equal to the number of final conditions.')
end
r = length(y0);			% Number of initial conditions
m = r + length(yf);		% Number of boundary conditions
% Checking the number of equations
ftest = feval(ODEfile,x0,[y0 ; guess],varargin{:});
if length(ftest) ~= m
    error(' The number of equations is not equal to the number of boundary conditions.')
end
fprintf('\n Integrating. Please wait.\n\n')
% Coefficients of the Legendre polynomial
for k = 0 : n/2
    cl(2*k+1) = (-1)^k * gamma(2*n-2*k+1) / ...
        (2^n * gamma(k+1) * gamma(n-k+1) * gamma(n-2*k+1));
    if k < n/2
        cl(2*k+2) = 0;
    end
end
zl = roots(cl);				% Roots of the Legendre polynomial
z = [-1; sort(zl); 1];		% Collocation points (z)
x = (xf-x0)*z/2+(xf+x0)/2;	% Collocation points (x)
% Bulding the vector of starting values of the dependent variables
[p,q] = RK(ODEfile,x0,xf,(xf-x0)/20,[y0 ; guess],2,varargin{:});
for k = 1:m
    y(k,:) = spline(p,q(k,:),x');
end
y(r+1:m,end) = yf(1:m-r);
% Building the matrix A
Q(:,1) = ones(n+2,1);
C(:,1) = zeros(n+2,1);
for i = 1:n+1
    Q(:,i+1) = x.^i;
    C(:,i+1) = i*x.^(i-1);
end
A = C*inv(Q);
for k = 1:m
    k1 = (k-1)*(n+2)+1;
    k2 = k1 + n+1;
    Am(k1:k2,k1:k2) = A;		% Building the matrix Am
    Y(k1:k2) = y(k,:);		% Building the vector Y
end
Y = Y';							% Make it a column vector
Y1 = Y * 1.1;
iter = 0;
maxiter = 100;
F = zeros(m*(n+2),1);
Fa = zeros(m*(n+2),1);
dY = zeros(m*(n+2),1);
position = [];	% Collocation points excluding boundary conditions
for k = 1:m
    if k <= r
        position = [position, (k-1)*(n+2)+[2:n+2] ];
    else
        position = [position, (k-1)*(n+2)+[1:n+1] ];
    end
end
% Newton's method
while max(abs(Y1 - Y)) > tol & iter < maxiter
    iter = iter + 1;
    fprintf(' Iteration %3d\n',iter)
    Y1 = Y;
    % Building the vector F
    for k = 1:n+2
        F(k : n+2 : (m-1)*(n+2)+k) = ...
            feval(ODEfile,x(k),Y(k : n+2 : (m-1)*(n+2)+k),varargin{:});
    end
    fnk = Am * Y - F;
    
    % Set dY for derivation
    for k = 1:m*(n+1)
        if Y(position(k)) ~= 0
            dY(position(k)) = Y(position(k)) / 100;
        else
            dY(position(k)) = 0.01;
        end
    end
    
    % Calculation of the Jacobian matrix
    for k = 1:m
        for kk = 1:n+1
            a = Y;
            nc = (k-1)*(n+1)+kk;
            a(position(nc)) = Y(position(nc)) + dY(position(nc));
            for kkk = 1:n+2
                Fa(kkk : n+2 : (m-1)*(n+2)+kkk) = ...
                    feval(ODEfile,x(kkk),a(kkk:n+2:(m-1)*(n+2)+kkk),varargin{:});
            end
            fnka = Am * a - Fa;
            jacob(:,nc) = (fnka(position) - fnk(position))...
                / dY(position(nc));
        end
    end
    % Next approximation of the roots
    if det(jacob) == 0
        Y(position) = Y(position) + max([abs(dY(position)); 1.1*tol]);
    else
        Y(position) = Y(position) - rho * inv(jacob) * fnk(position);
    end
end
% Rearranging the y's
for k = 1:m
    k1 = (k-1)*(n+2)+1;
    k2 = k1 + n+1;
    y(k,:) = Y(k1:k2)';
end
x = x';
if iter >= maxiter
    disp('Warning : Maximum iterations reached.')
end

Answer 1

[Voss]

https://www.mathworks.com/matlabcentral/answers/96005-why-do-i-get-the-error-unrecognized-function-or-variable

[erfan]

Runge-Kutta method

[Voss]

Cool.

Answer 2

[Star Strider]

Apparently you have to supply the ‘RK’ function referenced in:

[p,q] = RK(ODEfile,x0,xf,(xf-x0)/20,[y0 ; guess],2,varargin{:});

.

[erfan]

Runge-Kutta method

[erfan]

how

[Star Strider]

I understand that.

You need to supply that function (and the ‘ODEfile’ function) in order to run your code.

[erfan]

how

[Star Strider]

Check the documentation for the code you are using. Those functions should be available somewhere in the same source.

The code appears to be from a fairly ancient textbook:

% (c) N. Mostoufi & A. Constantinides

% January 1, 1999

so it is quite likely that it may not work in more recent MATLAB releases witthout being updated.

It would probably be better to use either bvp4c or bvp5c with your differential equations, once you have coded them. Both of these were part of MATLAB prior to R2006a, so there is probably no reason to use the code you posted.