function [total_In,total_theta,total_u_h,total_v_h]=proof_Gray_Scott(m,N,omega,L,nu,scaleEigs)

Tot_Cont_Steps=8;
dir=-1;
step_size=.1;
split_mesh=.38;
gamma=0.15;
igamma=intval('0.15');
r_star=.003; 

t=zeros(1,m+1)';
t(1:floor(m/2))=linspace(0,split_mesh,floor(m/2))';
t(floor(m/2):m+1)=linspace(split_mesh,1,m+2-floor(m/2))';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The starting point of the continuation is the exact homoclinic solution %
% given by THEOREM A for gamma in (0,2/9) on page 3 of                    %
%                                                                         %
% HALE, PELETIER and TROY, Exact homoclinic and heteroclinic solutions of %
% the Gray-Scott model for autocatalysis. SIAM Journal of Applied Math.,  %
% Vol. 61, No 1, pp 102-130, 2000.                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Q=sqrt(1-9*gamma/2);
u_h=1-3*gamma./(1+Q*cosh(L*t/sqrt(gamma)));
v_h=3./(1+Q*cosh(L*t/sqrt(gamma)));

%%% The following starting initial value of theta is somehow arbitrary
theta=[.4;-.06];

%%% We initialize the variable for the beginning of the continuation steps

Cont_Step=1;
total_u_h=zeros(m+1,1);
total_v_h=zeros(m+1,1);
total_theta=zeros(2,1);
total_In=zeros(2,1);

%%% Beginning of the continuation steps

while Cont_Step<=Tot_Cont_Steps
    
lambda=1/gamma+Cont_Step*dir*step_size;
ilambda=intval(lambda);
ip=intval([1 0 0 0]);
iscaleEigs=intval(scaleEigs);
iL=intval(L);
inu=intval(nu);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                  %%%
%%% PART I : PARAMETERIZATION OF THE STABLE MANIFOLD %%%
%%%            &  VERIFICATION OF THEOREM 4.2        %%%
%%%                                                  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Rigorous computation of the coefficients of the local stable manifold

%display('Rigorous Parameterization Evaluation');
[iA1,iA2,iA3,iA4,lambda_s1,lambda_s2]=intval_sym_grayScott_paramWs(ip,iscaleEigs,N,ilambda,igamma,iL);

%display('Aposteriori Validation of the Parameterization');
[goodComputation,delta]=intval_JP_aposteriori_numBound(N,iL,igamma,ilambda,inu,iA1,iA2,iA3,iA4,lambda_s1,lambda_s2);

if goodComputation==0
    display('The aposteriori validation of the parametrization failed !')
    return
end

A1=mid(iA1);
A2=mid(iA2);
A3=mid(iA3);
A4=mid(iA4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                            %%%
%%% PART II : NEWTON METHOD ON THE PROJECTION  %%%
%%%                                            %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Quantities necessary to start running Newton's method

% We evaluate F^{(m,N)}(theta,u_h,v_h)
f=F_GrayScott(theta,u_h,v_h,L,gamma,lambda,A1,A2,A3,A4,N,t);

% We evaluate the Jacobian matrix of at (theta,u_h,v_h):
% DF^{(m,N)}(theta,u_h,v_h).
A_m=DF_GrayScott(theta,u_h,v_h,L,gamma,lambda,A1,A2,A3,A4,N,t);

% We evaluate a numerical inverse of the Jacobian matrix DF^{(m,N)}(theta,u_h,v_h):
dagA_m=inv(A_m);

% We initialize Newton's method
k=0;
x1=[theta;u_h;v_h];
tol=5e-14;

%display('Newton iterations');
while (k<=15) && (norm(f)> tol),
    x1=x1-dagA_m*f;
    theta1=x1(1:2);
    u_h1=x1(3:m+3);
    v_h1=x1(m+4:2*m+4);
    f=F_GrayScott(theta1,u_h1,v_h1,L,gamma,lambda,A1,A2,A3,A4,N,t);
    A_m=DF_GrayScott(theta1,u_h1,v_h1,L,gamma,lambda,A1,A2,A3,A4,N,t);
    dagA_m=inv(A_m);
    k=k+1;                            %k: number of Newton's iterations
end

%%%% We store the values of the approximate solution (theta,u_h,v_h)
theta=x1(1:2);
u_h=x1(3:m+3);
v_h=x1(m+4:2*m+4);

%%%% We define interval values

i0=intval(0);
i1=intval(1);
i2=intval(2);
i3=intval(3);
i8=intval(8);
iN=intval(N);
itheta=intval(theta);
iu_h=intval(u_h);
iv_h=intval(v_h);
it=intval(t);
idelta=intval(delta);
iomega=intval(omega);
idagA_m=intval(dagA_m);

%display('Evaluation of the Jacobian using interval arithmetic');
iA_m=intval_DF_GrayScott(itheta,iu_h,iv_h,iL,igamma,ilambda,iA1,iA2,iA3,iA4,N,it);

if norm(mag(intval(eye(2*m+4))-iA_m*idagA_m))>=1
    display('The matrix dagA_m is not injective)');
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                  %%%
%%% PART III : CONSTRUCTION OF THE RADII POLYNOMIALS %%%
%%%            &  VERIFICATION OF THEOREM 2.6        %%%
%%%                                                  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            %
% Computation of the bound Y %
%                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%
%%%% Y_F %%%%
%%%%%%%%%%%%%

int_f=intval_F_GrayScott(itheta,iu_h,iv_h,iL,igamma,ilambda,iA1,iA2,iA3,iA4,N,it);
Y_F=mag(idagA_m*int_f)+mag(idelta*intval(abs(dagA_m))*intval(ones(2*m+4,1)));

%%%%%%%%%%%%%%%%%
%%%% Y_infty %%%%
%%%%%%%%%%%%%%%%%

%%% The following computes 
%%% ||u_h*v_h^2-lambda*(1-u_h)||_infty AND ||v_h-u_h*v_h^2||_infty 

s1=zeros(1,m);
s2=zeros(1,m);

for k=0:m-1
    u_k=infsup(min(u_h(k+1),u_h(k+2)),max(u_h(k+1),u_h(k+2)));
    v_k=infsup(min(v_h(k+1),v_h(k+2)),max(v_h(k+1),v_h(k+2)));
    s1(k+1)=mag(u_k*v_k^2-ilambda*(i1-u_k));
    s2(k+1)=mag(v_k-u_k*v_k^2);
end

Y_infty=mag((iL^2/i8)*intval(max(max(mag((it(2:end)-it(1:end-1)).^2.*intval(s1)')),max(mag((i1/igamma)*(it(2:end)-it(1:end-1)).^2.*intval(s2)')))));

%%%%%%%%%%%
%%%% Y %%%%
%%%%%%%%%%%

Y=[Y_F;Y_infty];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            %
% Computation of the bound Z %
%                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%
%%%% Z_F %%%%
%%%%%%%%%%%%%

sigma1=infsup(theta(1)-r_star,theta(1)+r_star);
sigma2=infsup(theta(2)-r_star,theta(2)+r_star);

hat_sigma1=intval(mag(sigma1/inu));
hat_sigma2=intval(mag(sigma2/inu));

if or(sup(hat_sigma1)>=1,sup(hat_sigma2)>=1),
    display('FAILURE!!! |hat_sigma1|>=1 or |hat_sigma2|>=1')
    return
end

%%% Computation of the LAMBDA^{(i)} (i=1,2)

powers_sum1=((2:N).*(1:N-1).*sigma1.^(0:N-2))'*(sigma2.^(0:N));
SUM(1,1)=intval(mag(sum(sum(powers_sum1.*iA1(3:end,:)))));
SUM(2,1)=intval(mag(sum(sum(powers_sum1.*iA2(3:end,:)))));
SUM(3,1)=intval(mag(sum(sum(powers_sum1.*iA3(3:end,:)))));
SUM(4,1)=intval(mag(sum(sum(powers_sum1.*iA4(3:end,:)))));

powers_sum2=((1:N).*sigma1.^(0:N-1))'*((1:N).*sigma2.^(0:N-1));
SUM(1,2)=intval(mag(sum(sum(powers_sum2.*iA1(2:end,2:end)))));
SUM(2,2)=intval(mag(sum(sum(powers_sum2.*iA2(2:end,2:end)))));
SUM(3,2)=intval(mag(sum(sum(powers_sum2.*iA3(2:end,2:end)))));
SUM(4,2)=intval(mag(sum(sum(powers_sum2.*iA4(2:end,2:end)))));

powers_sum3=(sigma1.^(0:N))'*((2:N).*(1:N-1).*sigma2.^(0:N-2));
SUM(1,3)=intval(mag(sum(sum(powers_sum3.*iA1(:,3:end)))));
SUM(2,3)=intval(mag(sum(sum(powers_sum3.*iA2(:,3:end)))));
SUM(3,3)=intval(mag(sum(sum(powers_sum3.*iA3(:,3:end)))));
SUM(4,3)=intval(mag(sum(sum(powers_sum3.*iA4(:,3:end)))));

infinite_sum1_term1=(hat_sigma2^(N+1)/(inu^2*hat_sigma1^2*(i1-hat_sigma2)))*sum(intval(2:N).*intval(1:N-1).*(hat_sigma1/hat_sigma2).^(2:N));
infinite_sum1_term2=((iN^2-iN)*hat_sigma1^2+(i2-i2*iN^2)*hat_sigma1+iN^2+iN)*hat_sigma1^(N-1)/(inu^2*(i1-hat_sigma1)^3*(i1-hat_sigma2));
infinite_sum1=infinite_sum1_term1+infinite_sum1_term2;

infinite_sum2_term1=(hat_sigma2^N/(inu^2*hat_sigma1*(i1-hat_sigma2)^2))*sum(intval(1:N-1).*((iN-intval(1:N-1))*(i1-hat_sigma2)+i1).*(hat_sigma1/hat_sigma2).^(1:N-1));
infinite_sum2_term2=(iN+hat_sigma1-iN*hat_sigma1)*hat_sigma1^(N-1)/(inu^2*(i1-hat_sigma1)^2*(i1-hat_sigma2)^2);
infinite_sum2=infinite_sum2_term1+infinite_sum2_term2;

infinite_sum3_term1=(hat_sigma1^(N+1)/(inu^2*hat_sigma2^2*(i1-hat_sigma1)))*sum(intval(2:N).*intval(1:N-1).*(hat_sigma2/hat_sigma1).^(2:N));
infinite_sum3_term2=((iN^2-iN)*hat_sigma2^2+(i2-i2*N^2)*hat_sigma2+iN^2+iN)*hat_sigma2^(N-1)/(inu^2*(i1-hat_sigma2)^3*(i1-hat_sigma1));
infinite_sum3=infinite_sum3_term1+infinite_sum3_term2;

infinite_sum=(infinite_sum1+i2*infinite_sum2+infinite_sum3)*idelta;

LAMBDA_0=[ SUM(1,1)+i2*SUM(1,2)+SUM(1,3)+infinite_sum ; SUM(3,1)+i2*SUM(3,2)+SUM(3,3)+infinite_sum];
LAMBDA_1=[ SUM(2,1)+i2*SUM(2,2)+SUM(2,3)+infinite_sum ; SUM(3,1)+i2*SUM(3,2)+SUM(3,3)+infinite_sum];

%%% We now use interval arithmetic to find upper bound for the Gamma^(ell)
%%% and tilde{Gamma}^(ell) for ell \in {1,2,3}
%%% (1) int_{t_k}^{t_{k+1}} |v^(l)(s)| ds, for l=1,2,3 and for 
%%% (2) int_{t_k}^{t_{k+1}} |(s-1) v^(l)(s)| ds, for l=1,2,3

sum_v1=intval(zeros(2,m));
sum_s_1v1=intval(zeros(2,m));
sum_v2=intval(zeros(2,m));
sum_s_1v2=intval(zeros(2,m));
sum_v3=intval(zeros(2,m));
sum_s_1v3=intval(zeros(2,m));

for k=0:m-1
    u_k=infsup(min(u_h(k+1),u_h(k+2)),max(u_h(k+1),u_h(k+2)));
    v_k=infsup(min(v_h(k+1),v_h(k+2)),max(v_h(k+1),v_h(k+2)));
    s=infsup(t(k+1),t(k+2));
    sum_v1(1,k+1)=(abs(v_k^2+ilambda)+abs(i2*u_k*v_k))*(it(k+2)-it(k+1))*iomega;
    sum_v1(2,k+1)=(i1/igamma)*(abs(v_k^2)+abs(i1-i2*u_k*v_k))*(it(k+2)-it(k+1))*iomega;
    sum_s_1v1(1,k+1)=(abs((s-i1)*(v_k^2+ilambda))+abs(i2*(s-i1)*u_k*v_k))*(it(k+2)-it(k+1))*iomega;
    sum_s_1v1(2,k+1)=(i1/igamma)*(abs((s-i1)*v_k^2)+abs((s-i1)*(i1-i2*u_k*v_k)))*(it(k+2)-it(k+1))*iomega;
    sum_v2(1,k+1)=i2*(abs(u_k)+abs(i2*v_k))*(it(k+2)-it(k+1))*(iomega+i1)^2;
    sum_v2(2,k+1)=(i2/igamma)*(abs(u_k)+abs(i2*v_k))*(it(k+2)-it(k+1))*(iomega+i1)^2;
    sum_s_1v2(1,k+1)=i2*(abs((s-i1)*u_k)+abs(i2*(s-i1)*v_k))*(it(k+2)-it(k+1))*(iomega+i1)^2;
    sum_s_1v2(2,k+1)=(i2/igamma)*(abs((s-i1)*u_k)+abs(i2*(s-i1)*v_k))*(it(k+2)-it(k+1))*(iomega+i1)^2;
    sum_v3(1,k+1)=i3*(it(k+2)-it(k+1))*(iomega+i1)^3;
    sum_v3(2,k+1)=(i3/igamma)*(it(k+2)-it(k+1))*(iomega+i1)^3;
    sum_s_1v3(1,k+1)=i3*((it(k+2)-it(k+1))*(i1-(it(k+1)+it(k+2))/i2))*(iomega+i1)^3;
    sum_s_1v3(2,k+1)=(i3/igamma)*((it(k+2)-it(k+1))*(i1-(it(k+1)+it(k+2))/i2))*(iomega+i1)^3;
end

%%%% Let us now construct the vectors V^(1), V^(2) and V^(3)

V1=intval(zeros(2*m+4,1));
V2=intval(zeros(2*m+4,1));
V3=intval(zeros(2*m+4,1));

V1(1,1)=iL^2*sum(sum_v1(1,:));
V1(2,1)=iL^2*sum(sum_v1(2,:));
V2(1,1)=LAMBDA_1(1)+iL^2*sum(sum_v2(1,:));
V2(2,1)=LAMBDA_1(2)+iL^2*sum(sum_v2(2,:));
V3(1,1)=iL^2*sum(sum_v3(1,:));
V3(2,1)=iL^2*sum(sum_v3(2,:));

V1(3,1)=iL^2*sum(sum_s_1v1(1,1:m-1));
V1(m+4,1)=iL^2*sum(sum_s_1v1(2,1:m-1));
V2(3,1)=LAMBDA_0(1)+iL^2*sum(sum_s_1v2(1,1:m-1));
V2(m+4,1)=LAMBDA_0(2)+iL^2*sum(sum_s_1v2(2,1:m-1));
V3(3,1)=iL^2*sum(sum_s_1v3(1,1:m-1));
V3(m+4,1)=iL^2*sum(sum_s_1v3(2,1:m-1));

for j=1:m-1
    V1(j+3,1)=(i1-it(j+1))*iL^2*sum(sum_v1(1,1:j))+iL^2*sum(sum_s_1v1(1,j:m-1));
    V1(m+j+4,1)=(i1-it(j+1))*iL^2*sum(sum_v1(2,1:j))+iL^2*sum(sum_s_1v1(2,j:m-1));
    V2(j+3,1)=LAMBDA_0(1)+(i1-it(j+1))*iL^2*sum(sum_v2(1,1:j))+iL^2*sum(sum_s_1v2(1,j:m-1));
    V2(m+j+4,1)=LAMBDA_0(2)+(i1-it(j+1))*iL^2*sum(sum_v2(2,1:j))+iL^2*sum(sum_s_1v2(2,j:m-1));
    V3(j+3,1)=(i1-it(j+1))*iL^2*sum(sum_v3(1,1:j))+iL^2*sum(sum_s_1v3(1,j:m-1));
    V3(m+j+4,1)=(i1-it(j+1))*iL^2*sum(sum_v3(2,1:j))+iL^2*sum(sum_s_1v3(2,j:m-1));
end

V1(m+3,1)=i0;
V1(2*m+4,1)=i0;
V2(m+3,1)=LAMBDA_0(1);
V2(2*m+4,1)=LAMBDA_0(2);
V3(m+3,1)=i0;
V3(2*m+4,1)=i0;

Z_F=zeros(2*m+4,3);
Z_F(:,1)=mag(intval(mag(intval(eye(2*m+4))-idagA_m*iA_m))*intval(ones(2*m+4,1)))+mag(intval(abs(dagA_m))*V1);
Z_F(:,2)=mag(intval(abs(dagA_m))*V2);
Z_F(:,3)=mag(intval(abs(dagA_m))*V3);

%%%%%%%%%%%%%%%%%
%%%% Z_infty %%%%
%%%%%%%%%%%%%%%%%

%%% The following computes 
%%% ||u_h*v_h^2-lambda*(1-u_h)||_infty AND ||v_h-u_h*v_h^2||_infty 

sZinfty1=zeros(1,m);
sZinfty2=zeros(1,m);
sZinfty3=zeros(1,m);

for k=0:m-1
    u_k=infsup(min(u_h(k+1),u_h(k+2)),max(u_h(k+1),u_h(k+2)));
    v_k=infsup(min(v_h(k+1),v_h(k+2)),max(v_h(k+1),v_h(k+2)));
    sZinfty1(k+1)=max(mag((abs(v_k^2+ilambda)+abs(i2*u_k*v_k))*(iomega+i1)),mag((i1/igamma)*(abs(v_k^2)+abs(i1-i2*u_k*v_k))*(iomega+i1)));
    sZinfty2(k+1)=max(mag(i2*(abs(u_k)+abs(i2*v_k))*(iomega+i1)^2),mag((i2/igamma)*(abs(u_k)+abs(i2*v_k))*(iomega+i1)^2));
    sZinfty3(k+1)=max(mag(i3*(iomega+i1)^3),mag((i3/igamma)*(iomega+i1)^3));
end

Z_infty(1,1)=mag(iL^2/i8*intval(max(mag((it(2:end)-it(1:end-1)).^2.*intval(sZinfty1')))));
Z_infty(1,2)=mag(iL^2/i8*intval(max(mag((it(2:end)-it(1:end-1)).^2.*intval(sZinfty2')))));
Z_infty(1,3)=mag(iL^2/i8*intval(max(mag((it(2:end)-it(1:end-1)).^2.*intval(sZinfty3')))));

%%%%%%%%%%%
%%%% Z %%%%
%%%%%%%%%%%

Z=[Z_F;Z_infty];

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Radii polynomials %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

ball_radii=[ones(1,2*m+4) omega];
p(1:2*m+5,1:4) = [Y Z];
p(1:2*m+5,2)=p(1:2*m+5,2)-ball_radii';
In=[0,Inf];
k=-1;

while In(2)>In(1) && k<=2*m+3,
  R=roots(p(k+2,4:-1:1));
  if norm(imag(R))>0
      display(['Failure with the ' ,num2str(k+2) 'th radii polynomial ' ]);
      return
  else
      if max(R)<0
          display(['Failure with the ' ,num2str(k+2) 'th radii polynomial ' ]);
          return
      end
  [~,n_min]=min(R);
  R(n_min)=[];
  I_k(1)=min(R);
  I_k(2)=max(R);
  In(1)=max(In(1),I_k(1));
  In(2)=min(In(2),I_k(2));
  k=k+1;
  end
end

if In(1)>=In(2)      
    display('Empty Intersection')
    return
elseif min(In)>r_star
    display('FAILURE!!! r_star > inf(In)')
    return
else
    %%% We really make sure that the polynomials are negative at the radius
    r0=(In(1)+In(2))/2;
    int_p(1:2*m+5,1:4) = [intval(Y) intval(Z)];
    int_p(1:2*m+5,2)=int_p(1:2*m+5,2)-intval(ball_radii');
    pr=int_p(:,4);
    for j=3:-1:1
        pr=int_p(:,j)+intval(r0)*pr;
    end
    if max(sup(pr))>=0
        display('FAILURE at proving that the polynomials are really negative')
        return
    end
    display(['Success with the validation at lambda = ' ,num2str(lambda)]) 
end

if max(mag(infsup(theta(1)-In(1),theta(1)+In(1))),mag(infsup(theta(2)-In(1),theta(2)+In(1))))>nu
    display('FAILURE!!! Theta is not in B_nu(0)')
    return
end

total_theta(:,Cont_Step)=theta;
total_u_h(:,Cont_Step)=u_h;
total_v_h(:,Cont_Step)=v_h;
total_In(:,Cont_Step)=In';

Cont_Step=Cont_Step+1;

end

end

function [F]=F_GrayScott(theta,u,v,L,gamma,lambda,A1,A2,A3,A4,N,t)

m=length(u)-1;
theta1=theta(1);
theta2=theta(2);

%%%%%%% We now evaluate the value of the parameterization of the manifolds

theta_powers=(theta1.^(0:N))'*(theta2.^(0:N));
P1_theta=sum(sum(theta_powers.*A1));
P2_theta=sum(sum(theta_powers.*A2));
P3_theta=sum(sum(theta_powers.*A3));
P4_theta=sum(sum(theta_powers.*A4));

F=zeros(2*m+4,1);

%%% Computation of the integrals on each of the mesh intervals [t_k,t_{k+1}] for k=0,...,m-1
I1=(1/12.*(t(2:m+1)-t(1:m))).*(6.*lambda.*u(1:m)+v(2:m+1).^2.*u(1:m)+2.*v(1:m).*v(2:m+1).*u(1:m)+3.*v(1:m).^2.*u(1:m)+3.*v(2:m+1).^2.*u(2:m+1)+2.*v(1:m).*v(2:m+1).*u(2:m+1)+6.*lambda.*u(2:m+1)+v(1:m).^2.*u(2:m+1)-12.*lambda);
I2=-(1/12).*(t(2:m+1)-t(1:m)).*(v(2:m+1).^2.*u(1:m)+2.*v(1:m).*v(2:m+1).*u(1:m)+3.*v(1:m).^2.*u(1:m)+v(1:m).^2.*u(2:m+1)+2.*v(1:m).*v(2:m+1).*u(2:m+1)+3.*v(2:m+1).^2.*u(2:m+1)-6.*v(2:m+1)-6.*v(1:m))/gamma;
I3=(1/60.*(t(2:m+1)-t(1:m))).*(-10.*v(1:m).*v(2:m+1).*u(1:m)-15.*v(1:m).^2.*u(1:m)-5.*v(2:m+1).^2.*u(1:m)-30.*lambda.*u(1:m)-10.*v(1:m).*v(2:m+1).*u(2:m+1)-15.*v(2:m+1).^2.*u(2:m+1)-5.*v(1:m).^2.*u(2:m+1)-30.*lambda.*u(2:m+1)+(6.*v(1:m).*v(2:m+1).*u(1:m)+2.*v(2:m+1).^2.*u(1:m)+20.*lambda.*u(1:m)+12.*v(1:m).^2.*u(1:m)+3.*v(1:m).^2.*u(2:m+1)+3.*v(2:m+1).^2.*u(2:m+1)+4.*v(1:m).*v(2:m+1).*u(2:m+1)+10.*lambda.*u(2:m+1)-30.*lambda).*t(1:m)+(10.*lambda.*u(1:m)+3.*v(2:m+1).^2.*u(1:m)+3.*v(1:m).^2.*u(1:m)+4.*v(1:m).*v(2:m+1).*u(1:m)+20.*lambda.*u(2:m+1)+2.*v(1:m).^2.*u(2:m+1)+12.*v(2:m+1).^2.*u(2:m+1)+6.*v(1:m).*v(2:m+1).*u(2:m+1)-30.*lambda).*t(2:m+1)+60.*lambda);
I4=-(1/60).*(t(2:m+1)-t(1:m)).*(-10.*v(1:m).*v(2:m+1).*u(1:m)-15.*v(1:m).^2.*u(1:m)-5.*v(2:m+1).^2.*u(1:m)-10.*v(1:m).*v(2:m+1).*u(2:m+1)-15.*v(2:m+1).^2.*u(2:m+1)-5.*v(1:m).^2.*u(2:m+1)+(3.*v(1:m).^2.*u(1:m)+4.*v(1:m).*v(2:m+1).*u(1:m)+3.*v(2:m+1).^2.*u(1:m)+6.*v(1:m).*v(2:m+1).*u(2:m+1)+2.*v(1:m).^2.*u(2:m+1)+12.*v(2:m+1).^2.*u(2:m+1)-20.*v(2:m+1)-10.*v(1:m)).*t(2:m+1)+30.*v(2:m+1)+(2.*v(2:m+1).^2.*u(1:m)+12.*v(1:m).^2.*u(1:m)+6.*v(1:m).*v(2:m+1).*u(1:m)+3.*v(2:m+1).^2.*u(2:m+1)+3.*v(1:m).^2.*u(2:m+1)+4.*v(1:m).*v(2:m+1).*u(2:m+1)-10.*v(2:m+1)-20.*v(1:m)).*t(1:m)+30.*v(1:m))/gamma;

%%% We build the vectors of partial sums to evaluate the integrals
ps_I1=zeros(m+1,1); ps_I2=zeros(m+1,1); ps_I3=zeros(m+1,1); ps_I4=zeros(m+1,1);

ps_I3(1)=sum(I3);
ps_I4(1)=sum(I4);

for j=1:m
	ps_I1(j+1)=ps_I1(j)+I1(j);
	ps_I2(j+1)=ps_I2(j)+I2(j);
	ps_I3(j+1)=ps_I3(j)-I3(j);
	ps_I4(j+1)=ps_I4(j)-I4(j);
end

F(1,1)=P2_theta-L^2*ps_I1(m+1);
F(2,1)=P4_theta-L^2*ps_I2(m+1);

F(3:m+3,1)=P1_theta+L^2*(t-1).*ps_I1+L^2*ps_I3-u;
F(m+4:2*m+4,1)=P3_theta+L^2*(t-1).*ps_I2+L^2*ps_I4-v;

end

function [DF]=DF_GrayScott(theta,u,v,L,gamma,lambda,A1,A2,A3,A4,N,t)

m=length(u)-1;
DF=zeros(2*m+4);

%%% Quantities required to compute the derivatives involving 
%%% the coefficients of the parameterization of the manifolds

theta_powers1=((0:N).*(theta(1).^[0 (0:N-1)]))'*(theta(2).^(0:N));
theta_powers2=(theta(1).^(0:N))'*((0:N).*(theta(2).^[0 (0:N-1)]));

%%% We define the vectors Il1, Il2, Il3, Il4, Il5, Il6, Il7, Il8 such that 
%%% Ilj=(0,I_1^{(l,j)},I_2^{(l,j)},..,I_m^{(l,j)}) \in R^(m+1) (j=1,2,3,4,5,6,7,8)
%%%% Note that since the nontrivial indices are k=1,...,m.

Il1=[0;(1/12.*(t(2:m+1)-t(1:m))).*(v(1:m).^2+3.*v(2:m+1).^2+6.*lambda+2.*v(1:m).*v(2:m+1))]';
Il2=[0;(1/6.*(v(2:m+1).*u(1:m)+v(1:m).*u(1:m)+3.*v(2:m+1).*u(2:m+1)+v(1:m).*u(2:m+1))).*(t(2:m+1)-t(1:m))]';
Il3=[0;-(1/12).*(2.*v(1:m).*v(2:m+1)+v(1:m).^2+3.*v(2:m+1).^2).*(t(2:m+1)-t(1:m))/gamma]';
Il4=[0;-(1/6).*(v(2:m+1).*u(1:m)+v(1:m).*u(1:m)+v(1:m).*u(2:m+1)+3.*v(2:m+1).*u(2:m+1)-3).*(t(2:m+1)-t(1:m))/gamma]';
Il5=[0;(1/60.*(t(2:m+1)-t(1:m))).*((10.*lambda+4.*v(1:m).*v(2:m+1)+3.*v(1:m).^2+3.*v(2:m+1).^2).*t(1:m)+(6.*v(1:m).*v(2:m+1)+2.*v(1:m).^2+12.*v(2:m+1).^2+20.*lambda).*t(2:m+1)-10.*v(1:m).*v(2:m+1)-15.*v(2:m+1).^2-30.*lambda-5.*v(1:m).^2)]';
Il6=[0;(1/30.*(t(2:m+1)-t(1:m))).*(-5.*v(2:m+1).*u(1:m)-5.*v(1:m).*u(1:m)-15.*v(2:m+1).*u(2:m+1)-5.*v(1:m).*u(2:m+1)+(2.*v(2:m+1).*u(1:m)+3.*v(1:m).*u(1:m)+2.*v(1:m).*u(2:m+1)+3.*v(2:m+1).*u(2:m+1)).*t(1:m)+(2.*v(1:m).*u(1:m)+3.*v(2:m+1).*u(1:m)+12.*v(2:m+1).*u(2:m+1)+3.*v(1:m).*u(2:m+1)).*t(2:m+1))]';
Il7=[0;(1/60).*(t(1:m)-t(2:m+1)).*((3.*v(2:m+1).^2+3.*v(1:m).^2+4.*v(1:m).*v(2:m+1)).*t(1:m)+(6.*v(1:m).*v(2:m+1)+2.*v(1:m).^2+12.*v(2:m+1).^2).*t(2:m+1)-15.*v(2:m+1).^2-5.*v(1:m).^2-10.*v(1:m).*v(2:m+1))/gamma]';
Il8=[0;-(1/30).*(t(2:m+1)-t(1:m)).*(-5.*v(2:m+1).*u(1:m)-5.*v(1:m).*u(1:m)-5.*v(1:m).*u(2:m+1)-15.*v(2:m+1).*u(2:m+1)+(2.*v(2:m+1).*u(1:m)+3.*v(1:m).*u(1:m)+2.*v(1:m).*u(2:m+1)+3.*v(2:m+1).*u(2:m+1)-5).*t(1:m)+(3.*v(2:m+1).*u(1:m)+2.*v(1:m).*u(1:m)+3.*v(1:m).*u(2:m+1)+12.*v(2:m+1).*u(2:m+1)-10).*t(2:m+1)+15)/gamma]';

%%% We define the vectors Ir1, Ir2, Ir3, Ir4, Ir5, Ir6, Ir7, Ir8 such that 
%%% Irj=(I_0^{(r,j)},I_1^{(r,j)},..,I_{m-1}^{(r,j)},0) \in R^(m+1) (j=1,2,3,4,5,6,7,8)
%%%% Note that since the nontrivial indices are k=0,...,m-1.

Ir1=[(1/12.*(t(2:m+1)-t(1:m))).*(3.*v(1:m).^2+2.*v(1:m).*v(2:m+1)+6.*lambda+v(2:m+1).^2);0]';
Ir2=[(1/6.*(3.*v(1:m).*u(1:m)+v(2:m+1).*u(1:m)+v(1:m).*u(2:m+1)+v(2:m+1).*u(2:m+1))).*(t(2:m+1)-t(1:m));0]';
Ir3=[-(1/12).*(v(2:m+1).^2+2.*v(1:m).*v(2:m+1)+3.*v(1:m).^2).*(t(2:m+1)-t(1:m))/gamma;0]';
Ir4=[-(1/6).*(3.*v(1:m).*u(1:m)+v(2:m+1).*u(1:m)+v(1:m).*u(2:m+1)+v(2:m+1).*u(2:m+1)-3).*(t(2:m+1)-t(1:m))/gamma;0]';
Ir5=[(1/60.*(t(2:m+1)-t(1:m))).*((6.*v(1:m).*v(2:m+1)+2.*v(2:m+1).^2+12.*v(1:m).^2+20.*lambda).*t(1:m)+(10.*lambda+3.*v(2:m+1).^2+4.*v(1:m).*v(2:m+1)+3.*v(1:m).^2).*t(2:m+1)-5.*v(2:m+1).^2-15.*v(1:m).^2-10.*v(1:m).*v(2:m+1)-30.*lambda);0]';
Ir6=[(1/30.*(t(2:m+1)-t(1:m))).*(-15.*v(1:m).*u(1:m)-5.*v(2:m+1).*u(1:m)-5.*v(2:m+1).*u(2:m+1)-5.*v(1:m).*u(2:m+1)+(3.*v(2:m+1).*u(1:m)+12.*v(1:m).*u(1:m)+2.*v(2:m+1).*u(2:m+1)+3.*v(1:m).*u(2:m+1)).*t(1:m)+(3.*v(1:m).*u(1:m)+2.*v(2:m+1).*u(1:m)+3.*v(2:m+1).*u(2:m+1)+2.*v(1:m).*u(2:m+1)).*t(2:m+1));0]';
Ir7=[(1/60).*(t(1:m)-t(2:m+1)).*((6.*v(1:m).*v(2:m+1)+2.*v(2:m+1).^2+12.*v(1:m).^2).*t(1:m)+(3.*v(1:m).^2+3.*v(2:m+1).^2+4.*v(1:m).*v(2:m+1)).*t(2:m+1)-10.*v(1:m).*v(2:m+1)-5.*v(2:m+1).^2-15.*v(1:m).^2)/gamma;0]';
Ir8=[-(1/30).*(t(2:m+1)-t(1:m)).*(-5.*v(2:m+1).*u(1:m)-15.*v(1:m).*u(1:m)-5.*v(1:m).*u(2:m+1)-5.*v(2:m+1).*u(2:m+1)+(3.*v(2:m+1).*u(1:m)+12.*v(1:m).*u(1:m)+2.*v(2:m+1).*u(2:m+1)+3.*v(1:m).*u(2:m+1)-10).*t(1:m)+(3.*v(1:m).*u(1:m)+2.*v(2:m+1).*u(1:m)+2.*v(1:m).*u(2:m+1)+3.*v(2:m+1).*u(2:m+1)-5).*t(2:m+1)+15)/gamma;0]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The first and second rows J_{1,*} and J_{2,*} %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DF(1,1)=sum(sum(theta_powers1.*A2));
DF(1,2)=sum(sum(theta_powers2.*A2));
DF(1,3:m+3)=-L^2*(Il1+Ir1);
DF(1,m+4:2*m+4)=-L^2*(Il2+Ir2);

DF(2,1)=sum(sum(theta_powers1.*A4));
DF(2,2)=sum(sum(theta_powers2.*A4));
DF(2,3:m+3)=-L^2*(Il3+Ir3);
DF(2,m+4:2*m+4)=-L^2*(Il4+Ir4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The first two columns of J %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DF(3:m+3,1)=sum(sum(theta_powers1.*A1));
DF(m+4:2*m+4,1)=sum(sum(theta_powers1.*A3));
DF(3:m+3,2)=sum(sum(theta_powers2.*A1));
DF(m+4:2*m+4,2)=sum(sum(theta_powers2.*A3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% We built the four sub-matrices %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=0:m

%%% Computation of J(3:m+3,3:m+3)
DF(3:k+2,k+3)=L^2*(Il5(k+1)+Ir5(k+1));
DF(k+3,k+3)=L^2*(t(k+1)-1)*Il1(k+1)+L^2*Ir5(k+1)-1;
DF(k+4:m+3,k+3)=L^2*(t(k+2:m+1)-1)*(Il1(k+1)+Ir1(k+1));

%%% Computation of J(3:m+3,m+4:2m+4)
DF(3:k+2,k+m+4)=L^2*(Il6(k+1)+Ir6(k+1));
DF(k+3,k+m+4)=L^2*(t(k+1)-1)*Il2(k+1)+L^2*Ir6(k+1);
DF(k+4:m+3,k+m+4)=L^2*(t(k+2:m+1)-1)*(Il2(k+1)+Ir2(k+1));

%%% Computation of J(m+4:2m+4,3:m+3)

DF(m+4:k+m+3,k+3)=L^2*(Il7(k+1)+Ir7(k+1));
DF(k+m+4,k+3)=L^2*(t(k+1)-1)*Il3(k+1)+L^2*Ir7(k+1);
DF(k+m+5:2*m+4,k+3)=L^2*(t(k+2:m+1)-1)*(Il3(k+1)+Ir3(k+1));

%%% Computation of J(m+4:2m+4,m+4:2m+4)
DF(m+4:k+m+3,k+m+4)=L^2*(Il8(k+1)+Ir8(k+1));
DF(k+m+4,k+m+4)=L^2*(t(k+1)-1)*Il4(k+1)+L^2*Ir8(k+1)-1;
DF(k+m+5:2*m+4,k+m+4)=L^2*(t(k+2:m+1)-1)*(Il4(k+1)+Ir4(k+1));

end

end

function [F]=intval_F_GrayScott(theta,u,v,L,gamma,lambda,A1,A2,A3,A4,N,t)

m=length(u)-1;
theta1=theta(1);
theta2=theta(2);
i1=intval(1); i2=intval(2); i3=intval(3); i4=intval(4);
i5=intval(5); i6=intval(6); i10=intval(10); i12=intval(12);
i15=intval(15); i20=intval(20); i30=intval(30); i60=intval(60);

%%%%%%% We now evaluate the value of the parameterization of the manifolds

theta_powers=(theta1.^(0:N))'*(theta2.^(0:N));
P1_theta=sum(sum(theta_powers.*A1));
P2_theta=sum(sum(theta_powers.*A2));
P3_theta=sum(sum(theta_powers.*A3));
P4_theta=sum(sum(theta_powers.*A4));

F=intval(zeros(2*m+4,1));

%%% Computation of the integrals on each of the mesh intervals [t_k,t_{k+1}] for k=0,...,m-1
I1=(i1/i12.*(t(2:m+1)-t(1:m))).*(i6.*lambda.*u(1:m)+v(2:m+1).^2.*u(1:m)+i2.*v(1:m).*v(2:m+1).*u(1:m)+i3.*v(1:m).^2.*u(1:m)+i3.*v(2:m+1).^2.*u(2:m+1)+i2.*v(1:m).*v(2:m+1).*u(2:m+1)+i6.*lambda.*u(2:m+1)+v(1:m).^2.*u(2:m+1)-i12.*lambda);
I2=-(i1/i12).*(t(2:m+1)-t(1:m)).*(v(2:m+1).^2.*u(1:m)+i2.*v(1:m).*v(2:m+1).*u(1:m)+i3.*v(1:m).^2.*u(1:m)+v(1:m).^2.*u(2:m+1)+i2.*v(1:m).*v(2:m+1).*u(2:m+1)+i3.*v(2:m+1).^2.*u(2:m+1)-i6.*v(2:m+1)-i6.*v(1:m))/gamma;
I3=(i1/i60.*(t(2:m+1)-t(1:m))).*(-i10.*v(1:m).*v(2:m+1).*u(1:m)-i15.*v(1:m).^2.*u(1:m)-i5.*v(2:m+1).^2.*u(1:m)-i30.*lambda.*u(1:m)-i10.*v(1:m).*v(2:m+1).*u(2:m+1)-i15.*v(2:m+1).^2.*u(2:m+1)-i5.*v(1:m).^2.*u(2:m+1)-i30.*lambda.*u(2:m+1)+(i6.*v(1:m).*v(2:m+1).*u(1:m)+2.*v(2:m+1).^2.*u(1:m)+i20.*lambda.*u(1:m)+i12.*v(1:m).^2.*u(1:m)+i3.*v(1:m).^2.*u(2:m+1)+i3.*v(2:m+1).^2.*u(2:m+1)+i4.*v(1:m).*v(2:m+1).*u(2:m+1)+i10.*lambda.*u(2:m+1)-i30.*lambda).*t(1:m)+(i10.*lambda.*u(1:m)+i3.*v(2:m+1).^2.*u(1:m)+i3.*v(1:m).^2.*u(1:m)+i4.*v(1:m).*v(2:m+1).*u(1:m)+i20.*lambda.*u(2:m+1)+i2.*v(1:m).^2.*u(2:m+1)+i12.*v(2:m+1).^2.*u(2:m+1)+i6.*v(1:m).*v(2:m+1).*u(2:m+1)-i30.*lambda).*t(2:m+1)+i60.*lambda);
I4=-(i1/i60).*(t(2:m+1)-t(1:m)).*(-i10.*v(1:m).*v(2:m+1).*u(1:m)-i15.*v(1:m).^2.*u(1:m)-i5.*v(2:m+1).^2.*u(1:m)-i10.*v(1:m).*v(2:m+1).*u(2:m+1)-i15.*v(2:m+1).^2.*u(2:m+1)-i5.*v(1:m).^2.*u(2:m+1)+(i3.*v(1:m).^2.*u(1:m)+i4.*v(1:m).*v(2:m+1).*u(1:m)+i3.*v(2:m+1).^2.*u(1:m)+i6.*v(1:m).*v(2:m+1).*u(2:m+1)+i2.*v(1:m).^2.*u(2:m+1)+i12.*v(2:m+1).^2.*u(2:m+1)-i20.*v(2:m+1)-i10.*v(1:m)).*t(2:m+1)+i30.*v(2:m+1)+(i2.*v(2:m+1).^2.*u(1:m)+i12.*v(1:m).^2.*u(1:m)+i6.*v(1:m).*v(2:m+1).*u(1:m)+i3.*v(2:m+1).^2.*u(2:m+1)+i3.*v(1:m).^2.*u(2:m+1)+i4.*v(1:m).*v(2:m+1).*u(2:m+1)-i10.*v(2:m+1)-i20.*v(1:m)).*t(1:m)+i30.*v(1:m))/gamma;

%%% We build the vectors of partial sums to evaluate the integrals
ps_I1=intval(zeros(m+1,1)); ps_I2=intval(zeros(m+1,1)); 
ps_I3=intval(zeros(m+1,1)); ps_I4=intval(zeros(m+1,1));

ps_I3(1)=sum(I3);
ps_I4(1)=sum(I4);

for j=1:m
	ps_I1(j+1)=ps_I1(j)+I1(j);
	ps_I2(j+1)=ps_I2(j)+I2(j);
	ps_I3(j+1)=ps_I3(j)-I3(j);
	ps_I4(j+1)=ps_I4(j)-I4(j);
end

F(1,1)=P2_theta-L^2*ps_I1(m+1);
F(2,1)=P4_theta-L^2*ps_I2(m+1);

F(3:m+3,1)=P1_theta+L^2*(t-1).*ps_I1+L^2*ps_I3-u;
F(m+4:2*m+4,1)=P3_theta+L^2*(t-1).*ps_I2+L^2*ps_I4-v;

end

function [DF]=intval_DF_GrayScott(theta,u,v,L,gamma,lambda,A1,A2,A3,A4,N,t)

m=length(u)-1;
DF=intval(zeros(2*m+4));

%%% Quantities required to compute the derivatives involving 
%%% the coefficients of the parameterization of the manifolds

theta_powers1=(intval(0:N).*(theta(1).^[0 (0:N-1)]))'*(theta(2).^(0:N));
theta_powers2=(theta(1).^(0:N))'*(intval(0:N).*(theta(2).^[0 (0:N-1)]));

i0=intval(0);
i1=intval(1);
i2=intval(2);
i3=intval(3);
i4=intval(4);
i5=intval(5);
i6=intval(6);
i10=intval(10);
i12=intval(12);
i15=intval(15);
i20=intval(20);
i30=intval(30);
i60=intval(60);

%%% We define the vectors Il1, Il2, Il3, Il4, Il5, Il6, Il7, Il8 such that 
%%% Ilj=(0,I_1^{(l,j)},I_2^{(l,j)},..,I_m^{(l,j)}) \in R^(m+1) (j=1,2,3,4,5,6,7,8)
%%%% Note that since the nontrivial indices are k=1,...,m.

Il1=[i0;(i1/i12.*(t(2:m+1)-t(1:m))).*(v(1:m).^2+i3.*v(2:m+1).^2+i6.*lambda+i2.*v(1:m).*v(2:m+1))]';
Il2=[i0;(i1/i6.*(v(2:m+1).*u(1:m)+v(1:m).*u(1:m)+i3.*v(2:m+1).*u(2:m+1)+v(1:m).*u(2:m+1))).*(t(2:m+1)-t(1:m))]';
Il3=[i0;-(i1/i12).*(i2.*v(1:m).*v(2:m+1)+v(1:m).^2+i3.*v(2:m+1).^2).*(t(2:m+1)-t(1:m))/gamma]';
Il4=[i0;-(i1/i6).*(v(2:m+1).*u(1:m)+v(1:m).*u(1:m)+v(1:m).*u(2:m+1)+i3.*v(2:m+1).*u(2:m+1)-i3).*(t(2:m+1)-t(1:m))/gamma]';
Il5=[i0;(i1/i60.*(t(2:m+1)-t(1:m))).*((i10.*lambda+i4.*v(1:m).*v(2:m+1)+i3.*v(1:m).^2+i3.*v(2:m+1).^2).*t(1:m)+(i6.*v(1:m).*v(2:m+1)+i2.*v(1:m).^2+i12.*v(2:m+1).^2+i20.*lambda).*t(2:m+1)-i10.*v(1:m).*v(2:m+1)-i15.*v(2:m+1).^2-i30.*lambda-i5.*v(1:m).^2)]';
Il6=[i0;(i1/i30.*(t(2:m+1)-t(1:m))).*(-i5.*v(2:m+1).*u(1:m)-i5.*v(1:m).*u(1:m)-i15.*v(2:m+1).*u(2:m+1)-i5.*v(1:m).*u(2:m+1)+(i2.*v(2:m+1).*u(1:m)+i3.*v(1:m).*u(1:m)+i2.*v(1:m).*u(2:m+1)+i3.*v(2:m+1).*u(2:m+1)).*t(1:m)+(i2.*v(1:m).*u(1:m)+i3.*v(2:m+1).*u(1:m)+i12.*v(2:m+1).*u(2:m+1)+i3.*v(1:m).*u(2:m+1)).*t(2:m+1))]';
Il7=[i0;(i1/i60).*(t(1:m)-t(2:m+1)).*((i3.*v(2:m+1).^2+i3.*v(1:m).^2+i4.*v(1:m).*v(2:m+1)).*t(1:m)+(i6.*v(1:m).*v(2:m+1)+i2.*v(1:m).^2+i12.*v(2:m+1).^2).*t(2:m+1)-i15.*v(2:m+1).^2-i5.*v(1:m).^2-i10.*v(1:m).*v(2:m+1))/gamma]';
Il8=[i0;-(i1/i30).*(t(2:m+1)-t(1:m)).*(-i5.*v(2:m+1).*u(1:m)-i5.*v(1:m).*u(1:m)-i5.*v(1:m).*u(2:m+1)-i15.*v(2:m+1).*u(2:m+1)+(i2.*v(2:m+1).*u(1:m)+i3.*v(1:m).*u(1:m)+i2.*v(1:m).*u(2:m+1)+i3.*v(2:m+1).*u(2:m+1)-i5).*t(1:m)+(i3.*v(2:m+1).*u(1:m)+i2.*v(1:m).*u(1:m)+i3.*v(1:m).*u(2:m+1)+i12.*v(2:m+1).*u(2:m+1)-i10).*t(2:m+1)+i15)/gamma]';

%%% We define the vectors Ir1, Ir2, Ir3, Ir4, Ir5, Ir6, Ir7, Ir8 such that 
%%% Irj=(I_0^{(r,j)},I_1^{(r,j)},..,I_{m-1}^{(r,j)},0) \in R^(m+1) (j=1,2,3,4,5,6,7,8)
%%%% Note that since the nontrivial indices are k=0,...,m-1.

Ir1=[(i1/i12.*(t(2:m+1)-t(1:m))).*(i3.*v(1:m).^2+i2.*v(1:m).*v(2:m+1)+i6.*lambda+v(2:m+1).^2);i0]';
Ir2=[(i1/i6.*(i3.*v(1:m).*u(1:m)+v(2:m+1).*u(1:m)+v(1:m).*u(2:m+1)+v(2:m+1).*u(2:m+1))).*(t(2:m+1)-t(1:m));i0]';
Ir3=[-(i1/i12).*(v(2:m+1).^2+i2.*v(1:m).*v(2:m+1)+i3.*v(1:m).^2).*(t(2:m+1)-t(1:m))/gamma;i0]';
Ir4=[-(i1/i6).*(i3.*v(1:m).*u(1:m)+v(2:m+1).*u(1:m)+v(1:m).*u(2:m+1)+v(2:m+1).*u(2:m+1)-i3).*(t(2:m+1)-t(1:m))/gamma;i0]';
Ir5=[(i1/i60.*(t(2:m+1)-t(1:m))).*((i6.*v(1:m).*v(2:m+1)+i2.*v(2:m+1).^2+i12.*v(1:m).^2+i20.*lambda).*t(1:m)+(i10.*lambda+i3.*v(2:m+1).^2+i4.*v(1:m).*v(2:m+1)+i3.*v(1:m).^2).*t(2:m+1)-i5.*v(2:m+1).^2-i15.*v(1:m).^2-i10.*v(1:m).*v(2:m+1)-i30.*lambda);i0]';
Ir6=[(i1/i30.*(t(2:m+1)-t(1:m))).*(-i15.*v(1:m).*u(1:m)-i5.*v(2:m+1).*u(1:m)-i5.*v(2:m+1).*u(2:m+1)-i5.*v(1:m).*u(2:m+1)+(i3.*v(2:m+1).*u(1:m)+i12.*v(1:m).*u(1:m)+i2.*v(2:m+1).*u(2:m+1)+i3.*v(1:m).*u(2:m+1)).*t(1:m)+(i3.*v(1:m).*u(1:m)+i2.*v(2:m+1).*u(1:m)+i3.*v(2:m+1).*u(2:m+1)+i2.*v(1:m).*u(2:m+1)).*t(2:m+1));i0]';
Ir7=[(i1/i60).*(t(1:m)-t(2:m+1)).*((i6.*v(1:m).*v(2:m+1)+i2.*v(2:m+1).^2+i12.*v(1:m).^2).*t(1:m)+(i3.*v(1:m).^2+i3.*v(2:m+1).^2+i4.*v(1:m).*v(2:m+1)).*t(2:m+1)-i10.*v(1:m).*v(2:m+1)-i5.*v(2:m+1).^2-i15.*v(1:m).^2)/gamma;i0]';
Ir8=[-(i1/i30).*(t(2:m+1)-t(1:m)).*(-i5.*v(2:m+1).*u(1:m)-i15.*v(1:m).*u(1:m)-i5.*v(1:m).*u(2:m+1)-i5.*v(2:m+1).*u(2:m+1)+(i3.*v(2:m+1).*u(1:m)+i12.*v(1:m).*u(1:m)+i2.*v(2:m+1).*u(2:m+1)+i3.*v(1:m).*u(2:m+1)-i10).*t(1:m)+(i3.*v(1:m).*u(1:m)+i2.*v(2:m+1).*u(1:m)+i2.*v(1:m).*u(2:m+1)+i3.*v(2:m+1).*u(2:m+1)-i5).*t(2:m+1)+i15)/gamma;i0]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The first and second rows J_{1,*} and J_{2,*} %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DF(1,1)=sum(sum(theta_powers1.*A2));
DF(1,2)=sum(sum(theta_powers2.*A2));
DF(1,3:m+3)=-L^2*(Il1+Ir1);
DF(1,m+4:2*m+4)=-L^2*(Il2+Ir2);

DF(2,1)=sum(sum(theta_powers1.*A4));
DF(2,2)=sum(sum(theta_powers2.*A4));
DF(2,3:m+3)=-L^2*(Il3+Ir3);
DF(2,m+4:2*m+4)=-L^2*(Il4+Ir4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The first two columns of J %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DF(3:m+3,1)=sum(sum(theta_powers1.*A1));
DF(m+4:2*m+4,1)=sum(sum(theta_powers1.*A3));
DF(3:m+3,2)=sum(sum(theta_powers2.*A1));
DF(m+4:2*m+4,2)=sum(sum(theta_powers2.*A3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% We built the four sub-matrices %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=0:m

%%% Computation of J(3:m+3,3:m+3)
DF(3:k+2,k+3)=L^2*(Il5(k+1)+Ir5(k+1));
DF(k+3,k+3)=L^2*(t(k+1)-1)*Il1(k+1)+L^2*Ir5(k+1)-1;
DF(k+4:m+3,k+3)=L^2*(t(k+2:m+1)-1)*(Il1(k+1)+Ir1(k+1));

%%% Computation of J(3:m+3,m+4:2m+4)
DF(3:k+2,k+m+4)=L^2*(Il6(k+1)+Ir6(k+1));
DF(k+3,k+m+4)=L^2*(t(k+1)-1)*Il2(k+1)+L^2*Ir6(k+1);
DF(k+4:m+3,k+m+4)=L^2*(t(k+2:m+1)-1)*(Il2(k+1)+Ir2(k+1));

%%% Computation of J(m+4:2m+4,3:m+3)

DF(m+4:k+m+3,k+3)=L^2*(Il7(k+1)+Ir7(k+1));
DF(k+m+4,k+3)=L^2*(t(k+1)-1)*Il3(k+1)+L^2*Ir7(k+1);
DF(k+m+5:2*m+4,k+3)=L^2*(t(k+2:m+1)-1)*(Il3(k+1)+Ir3(k+1));

%%% Computation of J(m+4:2m+4,m+4:2m+4)
DF(m+4:k+m+3,k+m+4)=L^2*(Il8(k+1)+Ir8(k+1));
DF(k+m+4,k+m+4)=L^2*(t(k+1)-1)*Il4(k+1)+L^2*Ir8(k+1)-1;
DF(k+m+5:2*m+4,k+m+4)=L^2*(t(k+2:m+1)-1)*(Il4(k+1)+Ir4(k+1));

end

end
