syms t;
syms p;
A = [3, 1; -13, -3];
B = [1, 3];
pI_A = p .* eye(2) - A;
inv_pI_A = (pI_A)^(-1);
exp_A = ilaplace(inv_pI_A);

t0 = 0;
t1 = sym(pi/8);

x0_1 = [-1; 2];
x0_2 = [2; 1];

x1 = [3; -4];

new_x0_1 = x0_1 - subs(exp_A, t, t0 - t1) * x1;
new_x0_2 = x0_2 - subs(exp_A, t, t0 - t1) * x1;