# Issue for Implementing the Runge Kutta CQ method

I would like to solve convolution equation numerically by the so-called convolution quadrature (CQ) (see the picture) based on Runge-Kutta CQ method. However, I didn’t arrive to get the expected solution as in the attached picture (TextBook; Integral Equation Methods for Evolutionary PDE, page 139) with matlab code. You can see my attempt with julia where all the steps are true,

``````using TaylorSeries, LinearAlgebra, Plots

function weights( Ks::Function , N , T )
A = [(88 - 7*√6)/360 (296-169* √6)/1800  (-2+3* √6) / 225
(296+169* √6) /1800 (88+7*√6) /360  (-2-3* √6) / 225
(16-√6 )/36 (16+√6)/36  1/9 ] ;
b = [(16-√6)/36 , (16+√6)/36  , 1/9];
c = [(4-√6)/10 , (4 + √6)/10 , 1] ;

Δt = T/N;
m = length(c) ;
ω=Matrix{Vector{Float64}}(undef,m,m)

z   = Taylor1(Float64, N+1)

Δ(z) = inv( A + ( z / (1-z) ) * ones(m) * b' );

y = Ks.(Δ(z) / Δt);

for i in 1:3
for j in 1:3
ω[i,j] = y[i,j].coeffs
end
end

return ω, A,b,c;

end

function evalRKCQ( g::Function , Ks::Function , N , T)

Δt = T/N ; ts = (0:N) * Δt ;

ω ,  A , b , c = weights( Ks , N , T )

m = length(c)

gs = zeros(m , N+1) ;

for  j in 1:m
gs[j,:] =  g.(ts .+ Δt * c[j]);
end

U=zeros(m,N+1)

for i in 1:m
for j in 1:m

for k in 1:N+1

U[i,k] = sum( ω[i,j][k+1-n] * gs[j,n] for n in 1:k)

end
end
end

return U

end

N = 100; T = 4; ts = (0:N)*T/N;

Ks(s) = 1 /(1 - exp(-s));

g(t) = exp(-100*(t-0.5) ^2);

exact= @. g(ts)+g(ts-1)+g(ts-2)+g(ts-3)

U = evalRKCQ( g , Ks  , N , T)

plot(ts[2:N],exact[2:N])

scatter!(ts[2:N], U[3,1:N-1])
``````

This is my plot’s attempt