matlab symbolic integration with real gives a complex answer

I wish to get analytical (closed form) solution to the following integral in Matlab. However, Matlab gives me an answer with real & imaginary parts. How do I get it to produce an answer with just the "real" part. Here is the full code.

close all;
clear all;
clc;

syms t real;
syms thetak real;
syms sik real;
syms tbar real;
syms sjk real;

expr = exp(-thetak*((t-sik)^2 + (t-sjk)^2));
Bijk_raw = int(expr,t,0,1);
Bijk = simplify(collect(expand(Bijk_raw)));
fprintf('Bijk is as follows...\n');
pretty(Bijk);

Solution

The answer you obtain (if you have a Matlab version similar to mine) and that I reproduce here :

  /               /                     2 \ 
  |  1/2   1/2    |   thetak (sik - sjk)  | 
- | 2    pi    exp| - ------------------- | 
  \               \            2          / 
 /    /  1/2          1/2                 \ 
 |    | 2    (-thetak)    (sik i + sjk i) | 
 | erf| --------------------------------- | i - 
 \    \                 2                 / 

    /  1/2          1/2                       \   \ \ 
    | 2    (-thetak)    (sik i + sjk i - 2 i) |   | | 
 erf| --------------------------------------- | i | | / 
    \                    2                    /   / / 
             1/2 
 (4 (-thetak)   )

gives the impression that you have complex number i everywhere.

But in fact it is a false impression due to (-thetak)^(1/2).

Indeed, taking the square root of a negative number will generate a "i" which in turn will "kill" the other "i"s that are "in contact" with it. This cancellation will occur at different places due to the fact that (-thetak)^(1/2) can be found :

1) inside the erf expressions and

2) as a common denominator (last line).

Verify that rule i^2=-1 applies everywhere leaving no chance to the survival of any "i"...

Finally giving (I have set thetak=s^2 with s>0) :

  /               /                      \ 
  |  1/2   1/2    |   s^2 (sik - sjk)^2   |   
- | 2    pi    exp| - ------------------- | 
  \               \            2          / 
 /    /  1/2                   \ 
 |    | 2    s    (sik  + sjk ) | 
 | erf| ----------------------- |  - 
 \    \           2            / 

    /  1/2                          \   \ \ 
    | 2    s   (sik  + sjk  - 2 )   |   | | 
 erf| ----------------------------- |   | |   /   (4 s)
    \             2                 /   / /

Edit : You could have escaped integration. The idea is to convert the quadratic in $exp(-thetak*((t-sik)^2 + (t-sjk)^2))$ under the so-called "canonical form" which in your case is: $exp(-thetak*(((t-A)^2 + B))/C);$ where $A,B,C$ can be expressed as fonctions of sik and sjk (for example $A=(sik+sjk)/2$); in this way, setting $T=t-A$, you are brought back to a classical Gauss integral with formula :

$$\frac{2}/{\sqrt{\pi}} \int_a^b exp(-t^2} dt ) (erf(b) - erf(a))$$