Guessing about what you want:
function main
%%%%%%NUMERICAL CODE
A=0.5; pr=1; a=1; a1=2;a2=1;xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,10),[0 1 0 a 1 0 0 0]);
sol=bvp5c(@ode,@bc,solinit);
xn=linspace(xa,xb,100);
x = xn;
S=deval(sol,x);
function res=bc(ya,yb)
res=[ya(1); ya(2)-1; ya(4); ya(5)-a; ya(7)-1; yb(2); yb(5); yb(7)];
end
function dydx=ode(~,y)
dydx=[y(2); y(3); 2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));
y(5); y(6); 2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
y(8); A*pr*a2*y(7)*(y(2)+y(5))-pr*a2*y(8)*(y(1)+y(4))];
end
f0 = deval(sol,0);
p=f0(3);q=f0(6);r=f0(8);
figure(1)
plot(x,S([2],:)); %for f'
xlabel('\eta'); ylabel('f`');
hold on
% %%%%%%%%%%%%%%%%% SYMBOLIC CODE
t = sym('t');
x = sym('x');
%{
p = sym('p');
q = sym('q');
r = sym('r');
a1 = sym('a1');
a2 = sym('a2');
A = sym('A');
pr = sym('pr');
%}
f = zeros(1,3,'sym');
g = zeros(1,3,'sym');
h = zeros(1,3,'sym');
fa = zeros(1,3,'sym');
ga = zeros(1,3,'sym');
ha = zeros(1,3,'sym');
f(1)=x+(p/2)*x^2;g(1)=a*x+(q/2)*x^2;h(1)=1+r*x;
for i=1:3
fa(i) = subs(f(i),x,t);dfa = diff(fa(i),t,1);d2fa = diff(dfa,t,1);
ga(i) = subs(g(i),x,t);dga = diff(ga(i),t,1);d2ga = diff(dga,t,1);
ha(i) = subs(h(i),x,t);dha = diff(ha(i),t,1);
If1=int((fa(i)+ga(i))*(-d2fa) +2*dfa*(dfa+dga),t,0,t);If2=int(If1,t,0,t);If3=int(If2,t,0,x);
Ig1=int((fa(i)+ga(i))*(-d2ga) + 2*dfa*(dfa+dga),t,0,t);Ig2=int(Ig1,t,0,t);Ig3=int(Ig2,t,0,x);
Ih1=int((fa(i)+ga(i))*(-dha) + A*ha(i)*(dfa+dga),t,0,t);Ih2=int(Ih1,t,0,x);
f(i+1) = a1*If3;g(i+1) = a1*Ig3;h(i+1) = pr*a2*Ih2;
disp(f(i+1))
end
f=f(1)+f(2)+f(3);g=g(1)+g(2)+g(3);h=h(1)+h(2)+h(3);
F=[0 diff(double(subs(f,x,xn)))];
figure(2)
plot(xn,F,'LineWidth',1.5)
hold on
end
%%%%%%%%%
