BQN, 197 bytes

{a‿b‿c‿d:(b÷3×a)-˜•math{𝕩>0?+´𝕩(𝕗.Cbrt+⋈-)⟜√˜-q÷2;𝕩=0?0⊸=◶⟨¯1‿2‿2÷˜·𝕗.Cbrt×⟜4,3⊸⥊⟩q;(2×√-p÷3)×𝕗.Cos(2×π×↕⊸÷3)-˜3÷˜𝕗.Acos(√-3÷p)×1.5×q÷p}(27÷˜p⋆3)+4÷˜×˜q←(d÷a)-(27÷˜3⋆˜b÷a)+3÷˜b×a÷˜p←(c÷a)-3÷˜×˜b÷a}

BQN is a minimalistic language, it does not have built-ins to solve polynomial equations, nor does it have support for complex numbers (none of the available implementations, it is specified though). If a=0, the program returns NaN.

Explanation:

Cub ← {a‿b‿c‿d:                                         
  p ← (c÷a)-3÷˜×˜b÷a                                       # Coefficients of the depressed cubic equation
  q ← (d÷a)-(27÷˜3⋆˜b÷a)+3÷˜b×a÷˜p                         # t^3 + pt + q = 0
  r ← (27÷˜p⋆3)+4÷˜×˜q                                     # Discriminant, need to analyze three cases
  Solve ← •math{                                           # Use a 1-modifer to "wrap" the •math namespace   
    𝕩>0 ? +´𝕩(𝕗.Cbrt+⋈-)⟜√˜-q÷2;                           # Cardano's method
    𝕩=0 ? 0⊸=◶⟨¯1‿2‿2÷˜·𝕗.Cbrt×⟜4,3⊸⥊⟩q;                   # Cardano's method
    (2×√-p÷3)×𝕗.Cos(2×π×3÷˜↕3)-˜3÷˜𝕗.Acos(√-3÷p)×1.5×q÷p    # Viète's trigonometric method
  }
  (Solve r)-b÷3×a                                          # Convert back from depressed cubic
}

I have added test cases for the 4 branches, try it in the

BQN online REPL

APL(NARS), 202 chars

r←{1e¯8<∣11○⍵:3√⍵⋄⍵≥0:3√⍵⋄-3√-⍵}
f←{(a b c)←1↓⍵÷↑⍵⋄p←¯3÷⍨b-3÷⍨a*2⋄q←¯2÷⍨c+(27÷⍨2×a*3)-a×b÷3⋄(b c)←r¨q(+,-)√(q*2)-p*3⋄(-a÷3)++/¨(⊂b c),{⍵×e×1 e←2÷⍨¯1+0J1×√3}¨(c b)(b c)}
m←{{(∼⍵=∅)/⍵}{1e¯8>∣11○⍵:9○⍵⋄∅}¨f⍵}

Above function f calculate the one array of 3 roots of the equation of argument A B C D that means

Ax^3+Bx^2+Cx+D=0

m return only the real ones of that array and r make root3 of all Complex numbers in the way formula in function f works ok. In these functions I consider one number x ugual to 0 if abs(x)<1e¯8.

11○ and 9○ are the immaginary part and the real part of input.

If we have the equation

Ax^3+Bx^2+Cx+D=0

Because A!=0 we can divide for A and obtain

x^3+(B/A)x^2+(C/A)x+(D/A)=0

and so call (a b c) (B/A, C/A, D/A) we can rewrite the equation as

x^3+a x^2+b x+c=0    

If we have this change of variables x=y-(a/3) we can rewrite above equation as

y^3+y(b-(a^2)/3)+c-(ab)/3+(2/27)a^3=0

and from this is know a way of calculate solutions using radicals (from wikipedia) it will find y01 y02 y03 solution of above equation, and using the translation

x0i=y0i-(a/3)

we can find the roots of the equation x^3+a x^2+b x+c=0. This is an old and difficult problem in the history of mathematic connected with the theory of immaginary numbers.

From wikipedia on cubic + my math book it should be

x^3-3px-2q=0

if

b=3√(q+√(q^2-p^3))
c=3√(q-√(q^2-p^3))

where "3√" "means cubic root of" and "√" is the squareroot, has always the real solution (even if q^2-p^3<0)

x_1=b+c

and the 2 solution:

x_2=b*e+c*e^2
x_3=c*e+b*e^2

where e=(-1+i*√3)/2

Test:

  f 1 3 0 ¯1
0.5320888862J0 ¯0.6527036447J1.110223025E¯16 ¯2.879385242J¯5.551115123E¯17 
  m 1 3 0 ¯1
0.5320888862 ¯0.6527036447 ¯2.879385242 
  m 1 0 ¯1 ¯1
1.324717957 
  m 1 1 1 1
¯1 
  m 1 0 ¯3 ¯2
2 ¯1 ¯1 

MATLAB, 29 bytes

syms a b c d;roots([a b c d])

Wolfram Language (Mathematica), 35 bytes

#^0NSolve[x#+#2&~Fold~!##,x,Reals]&

Try it online!

Input [a, b, c, d]. If \$a=0\$, returns a list of as many {Indeterminate}s as there would be solutions.

When given non-equation expressions, NSolve finds their roots.

          x#+#2&~Fold~!##           convert to a polynomial in x
   NSolve[               ,x,Reals]  get real solutions
#^0                                 invalidate solutions if a=0

J, 20 bytes

1(#~]=+)@{::[:p.]*{:

Attempt This Online!

Takes d c b a, and returns the array of real roots, or an empty array if a is 0.

1(#~]=+)@{::[:p.]*{:
                ]*{:    multiply the entire polynomial with `a`
                        which doesn't change the roots if `a` is nonzero
                        but gives 0 0 0 0 otherwise, which has zero roots according to `p.`
            [:p.    convert to multipler-roots form
1        {::        extract the roots
 (     )@      filter real roots:
  #~]=+        filter the list by "complex conjugate equals self"

J, 32 bytes

[:{:^:([:-.*./@(=+)){:@:>@p.^:{:

I am pretty bad at conditionals in J, but after some finagling, this was the best I could come up with. Expects coefficients as d,c,b,a

Attempt This Online!

[:{:^:([:-.*./@(=+)){:@:>@p.^:{:
                    {:@:>@p.^:{:  : u^:v executes u if v is true, else return x
                              {:  : is last item 0?
                    {:@:>@p.      : if non-zero
                          p.      : converts coefficient form to multiplier-roots form
                        >         : unbox the results
                    {:            : take the roots, not the multiplier
  {:^:([:-.*./@(=+))              : The check for only one real
      ([:-.*./@(=+))              : u(v(x)) where u is -. and v is *./@(=+)
               (=+)               : each x equal to its complex conjugate?
           *./                    : AND reduce result
         -.                       : negate result
  {:                              : if roots are not all reals, take the last value

extra:
u@v -> apply u monadically to v(x) for every application of v
u@:v -> apply u to the entire result of v(x)
[:u v -> strictly u(v(x))

Python 3, 99 bytes

import numpy
f=lambda a,b,c,d:[float(i)for i in numpy.roots([a,b,c,d])if round(float(i))!=0if a!=0]

Try it online!

Jelly, 8? 11 bytes

ṛ/ȧÆr,N+AƲƇ

A monadic Link that accepts a list of the coefficients ([d,c,b,a]) and yields a list of the (one or three) real roots, unless a is zero in which case it yields an empty list.

Try it online! Or see the test-suite.

...\$8\$ bytes if we can handle a being zero (in which case it works for all polynomials):

Ær,NċAƲƇ

Try it online!

How?

The heavy lifting is done by a built-in...

ṛ/ȧÆr,N+AƲƇ - Link: coefficients, P = [d,b,c,a]
 /          - reduce (P) by:
ṛ           -   right
               -> a
   Ær       - polynomial roots (of P)
  ȧ         - logical AND
               -> [d,c,b,a] or 0
          Ƈ - keep those for which (with 0, keep those of range(0)=[]):
         Ʋ  -   last four links as a monad - f(x):
      N     -     negate (x)
     ,      -     (P) paired with (that)
        A   -     absolute value (of x)
       ċ    -     count occurrences (of that in the pair)
                    -> Truthy (1 or 2*) if x is real, falsey (0) if not
                                  * when x=0

Vyxal P, 4 bytes

hß∆P

Try It Online!

Good old built-ins. This may not work on some interpreters until the SymPy ACE vulnerabilities are fixed. Input as a list [a,b,c,d].

Edit: Changed output to the input if leading coefficient is 0.

MATL, 12 bytes

1)?G&ZQ&Zj~)

Try it online!

How it works

       % Implicit input
1)     % Get the first element
?      % If non-zero
  G    %   Push input again
  &ZQ  %   Implicit input. Roots of polynomial. Gives a numeric vector
  &Zj  %   Push real and imaginary parts
  ~    %   Negate. Gives "true" if imaginary part is zero, or "false" otherwise
  )    %   Keep only real parts corresponding to "true".
       % Implicit end
       % Implicit display

Python, 143 134 133 bytes

f=lambda a,b,c,d,R=[-1,1j,1]:[R:=[(z:=R[2])-(((a*z+b)*z+c)*z+d)/(z-R[0])/(z-R[1])]+R for i in" "*99]and[x for x in R[:3]if x.imag==0]

Attempt This Online!

This is the first solution that doesn't use a builtin, and it suffers greatly for it. Instead, I'm using the Durant-Kerner method to find all approximations, then filtering for whether they're complex.

-9 bytes from @Steffan for various small changes.

Wolfram Language (Mathematica), 50 bytes

Try it online!

g[k_]/;k[[4]]!=0:=Solve[x^Range[0,3].k==0,x,Reals]

Input is the number of coefficients in reverse order {d, c, b, a}. This defines a function g only if the 4th argument is unequal 0. Otherwise this function remains undefined and returns the function call as output. Uses the built-in Solve function over the Reals domain.

Both functions below take a polynomial in x.

PARI/GP, 27 bytes

p->if(#p>3,polrootsreal(p))

Attempt This Online!

Returns 0 for quadratics and below.

PARI/GP, 101 bytes

p->k=[Pi,I,2];for(q=1,99,k-=[subst(p,x,k[i])/vecprod([k[i]-z|z<-k[^i]])|i<-[1..3]]);[r|r<-k,!imag(r)]

Attempt This Online!

Yet another implementation of the Durand–Kerner method.

R, 42 41 bytes

function(c,a=polyroot(c))a[!Im(a)]/!!c[4]

Try it online!

Input is vector of d,c,b,a. Outputs Inf if the cubic coefficient (a) is zero.

Using a[!Im(a)] to select only real solutions is very susceptible to floating-point rounding errors; the TIO header rounds values less than 0.0000000001 to zero to prevent this. Including rounding in the code costs 5 3 more bytes.

R, 111 bytes

function(e,z=1i^(1:3)){for(i in 1:99)z=z-e%*%rbind(z^3,z^2,z,1)/(z-(x=c(z,z))[2:4])/(z-x[3:5]);if(e)z[!Im(z)]}	

Try it online!

Non-builtin Durant-Kerner strategy copied from 97.100.97.109's answer (upvote that!). Using this approach to find all the roots in parallel works very nicely with R's vectorized operations.