Solving a system of equations symbolically

I need to solve a system of equations.
Everything works until the part where I want to solve for the coefficients which I have stored in the vector a:
Is there a way to solve for these coefficients?

My script is the following:

using SymPy
using Symbolics

@variables a₁ a₂ a₃ a₄ L x
a = [a₁,a₂,a₃,a₄]

phi = 0;
for i in axes(a,1)
        phi = phi + a[i]*x^(i-1)
end

phid = Symbolics.derivative(phi, x)

# Defining the equations
equ1 = [substitute(phi,  x => 0) => 1,     # phi(0)=1
             substitute(phi,  x => L) => 0,     # phi(L)=0
             substitute(phid, x => 0) => 0,    # dphi(0)=0
             substitute(phid, x => L) => 0]    # dphi(L)=0

# The part that does not work!!::
# Solve for a₁ a₂ a₃ a₄
# sol1 = solve(equ1, a)
# sol1.a₁
# sol1.a₂
# sol1.a₃
# sol1.a₄


# Next step
# Substitute coefficients into phi
phi_t1 = substitute(phi, a => [sol1.a₁,sol1.a₂,sol1.a₃,sol1.a₄])

Not with Symbolics.jl yet. Use SymPy.jl for this for now.

Thank you for your response!
How would you then solve it using SymPy.jl?

Not very different than what you have:

using SymPy
@syms a[1:4] L x
phi = zero(x)
for i ∈ axes(a,1)
    phi = phi + a[i] * x^(i-1)
end

phid = diff(phi, x)
equ1 = (phi(x => 0) ~ 1,
        phi(x => L) ~ 0,
        phid(x => 0) ~ 0,
        phid(x => L) ~ 0)

solve(equ1, a)