function u = sol1Dwave(c,a,b)
x=linspace(-10*a,10*a,1000);
t=[0,a/(2*c),a/c,2*a/c,3*a/c,4*a/c,5*a/c];
T=length(t);
L=length(x);
Q=zeros(T,L);
u=zeros(T,L);
%phi1=zeros(1,L);
%phi2=zeros(1,L);
for i=1:T
 cr1=x+c*t(i);
 cr2=x-c*t(i);
 phi1 = initpos(cr1,b,a);
 phi2 = initpos(cr2,b,a);
 for j=1:L
 Q(i,j) = quad(@initmom,x(j)-c*t(i),x(j)+c*t(i));
 end
 u(i,:)=(1/2)*(phi1+phi2)+(1/(2*c))*Q(i,:);
end

figure(1)

subplot(6,1,1)
plot(x,u(1,:))
xlim([-10*a 10*a])
ylim([0,5])
title('Plot at t=0')
subplot(6,1,2)
plot(x,u(2,:))
xlim([-10*a 10*a])
ylim([0,5])
title('Plot at t=a/2c')
subplot(6,1,3)
plot(x,u(3,:))
xlim([-10*a 10*a])
ylim([0,5]);
title('Plot at t=a/c')
subplot(6,1,4)
plot(x,u(4,:))
xlim([-10*a 10*a])
ylim([0,5]);
title('Plot at t=2a/c')
subplot(6,1,5)
plot(x,u(5,:))
xlim([-10*a 10*a])
ylim([0,5]);
title('Plot at t=3a/c')
subplot(6,1,6)
plot(x,u(6,:))
xlim([-10*a 10*a])
ylim([0,5]);
title('Plot at t=4a/c')

figure(2)
surf(x,t,u)
 