function plotSolution(sol,y,t)
l=length(sol(:,1,1));
u1 = sol(l,:,1);
u2 = sol(l,:,2);
u3 = sol(l,:,3);
u4 = sol(l,:,4);
u5 = sol(l,:,5);
u6 = sol(l,:,6);
u7 = sol(l,:,7);
u8 = sol(l,:,8);
u9 = sol(l,:,9);
u10 = sol(l,:,10);
u11 = sol(l,:,11);
u12 = sol(l,:,12);
u13 = sol(l,:,13);
u14 = sol(l,:,14);
eS=sqrt((1+exp(-u14)).^(-1));
eS2=sqrt(u10.*u10+u11.*u11+u12.*u12);
alphaS=0.74.*(1+exp(-u13)).^(-1);
subplot(3,5,1)
s1=plot(u1,y);
title('Fluid velocity u_x(y,t)')
ylabel('y')
xlabel('u_x')
%s1.EdgeColor= 'none'

subplot(3,5,2)
s2=plot(u2,y);
title('Fluid velocity u_y(y,t)')
ylabel('y')
xlabel('u_y')

subplot(3,5,3)
s3=plot(u3,y);
title('Fluid velocity u_z(y,t)')
ylabel('y')
xlabel('u_z')

subplot(3,5,4)
s4=plot(u4,y);
title('Solid translational velocity v_x(y,t)')
ylabel('y')
xlabel('v_x')

subplot(3,5,5)
s5=plot(u5,y);
title('Solid translational velocity v_y(y,t)')
ylabel('y')
xlabel('v_y')

subplot(3,5,6)
s6=plot(u6,y);
title('Solid translational velocity v_z(y,t)')
ylabel('y')
xlabel('v_z')

subplot(3,5,7)
s7=plot(u7,y);
title('Solid rotational velocity omega_x(y,t)')
ylabel('y')
xlabel('omega_x')

subplot(3,5,8)
s8=plot(u8,y);
title('Solid rotational velocity omega_y(y,t)')
ylabel('y')
xlabel('omega_y')

subplot(3,5,9)
s9=plot(u9,y);
title('Solid rotational velocity omega_z(y,t)')
ylabel('y')
xlabel('omega_z')

subplot(3,5,10)
s10=plot(u10,y);
title('Particle orientation e_x(y,t)')
ylabel('y')
xlabel('e_x')

subplot(3,5,11)
s11=plot(u11,y);
title('Particle orientation e_y(y,t)')
ylabel('y')
xlabel('e_y')

subplot(3,5,12)
s12=plot(u12.*eS./eS2,y);
title('Particle orientation e_z(y,t)')
ylabel('y')
xlabel('e_z')

subplot(3,5,13)
s13=plot(eS,y);
title('Absolute value orientation vector |e|')
ylabel('y')
xlabel('e')

subplot(3,5,14)
s14=plot(alphaS.*(u1-u4),y);
title('Relative velocity w(y,t)=u_x-v_x')
ylabel('y')
xlabel('w')

subplot(3,5,15)
s14=plot(0.74.*(1+exp(-u13)).^(-1),y);
title('Solid phase volume fraction alphaS(y,t)')
ylabel('y')
xlabel('alphaS')

end