# Issue for Implementing the Runge Kutta CQ method

**URL:** <https://discourse.julialang.org/t/issue-for-implementing-the-runge-kutta-cq-method/97397>\
**Category:** General Usage\
**Tags:** question, plotting\
**Created:** [April 12, 2023, 3:27pm UTC](https://discourse.julialang.org/t/issue-for-implementing-the-runge-kutta-cq-method/97397 "2023-04-12T15:27:38Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [April 12, 2023, 3:27pm UTC](https://discourse.julialang.org/t/issue-for-implementing-the-runge-kutta-cq-method/97397/1 "2023-04-12T15:27:39Z")

</div>

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,

```julia
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])

```

 ![Capture2](https://global.discourse-cdn.com/julialang/original/3X/5/5/5594dd683ece2fc51cd8874d7b9aa67b17f4ea50.png)  
 ![Capture4](https://global.discourse-cdn.com/julialang/original/3X/3/d/3da4a80063e2205bdbca7420929d0851c89e9920.png)

This is my plot’s attempt

 ![Capture5](https://global.discourse-cdn.com/julialang/original/3X/a/6/a6d140ec656bcb36908bae05df107869c16543c5.png)