how to write the equation for expected prediction value of a gpr model in symbolic form

I am trying to get the prediction function of fitgpr in symbolic form.

My model type is Linear basis, exact FitMethod and exact PredictMethod. I kept the kernel function as default so it is squared exponential.

In this article, it says if using exact gpr method, i can use this equation for prediction

I coded(attached below) up the function in symbolic form and when i substitute xe to 0 to check the output, my value is off. I checked with other values of xe as well and my symbolic equation output came out wrong at all the values tested.

What am I doing wrong?

Here is my code

clear all
clc
close all

rng(0,'twister');   % For reproducibility

n = 11;               % n - size of dataset 
c = 5;                % c - drag co-efficient
xi = 0;               % xi - lower bound
xf = 5;               % xf - upper bound
var  = 0.2;           % v  - variance



x = linspace(xi,xf,n)';                     % populating x
y = (-c * x) + (var * randn(n,1));            % y - eqn

gprMdl = fitrgp(x,y,'Basis','Linear',...
    'FitMethod','exact','PredictMethod','exact');

syms v f

sig_l = gprMdl.KernelInformation.KernelParameters(1);
sig_f = gprMdl.KernelInformation.KernelParameters(2);

f = [v 1] * [gprMdl.Beta];

for i=1:n
    f = f + ((gprMdl.Alpha(i) * (sig_f^2)) * exp( -0.5 * (sig_l^(-2)) * ( v - x(i) ) ) ) ;
end

eval(subs(f,0))

predict(gprMdl,0)